\section{Introduction}
The present article is mostly motivated by the work \cite{Taborite} about negligible sets for convexity of functions in $\er^d$, where an interesting open problem was raised. We shall need the following notion introduced in \cite{Taborite}.

A set $A\subset \er^d$ is called \emph{intervally thin} if for any $x,y\in \er^d$ and any $\eps >0$ there exist $x'\in B(x,\eps)$ and $y'\in B(y,\eps)$ such that $[x',y']\cap A=\emptyset$.

\begin{problemTT}
Let $A\subset \er^n$ be closed. Suppose that for an arbitrary open set $U$ containing $A$ every locally convex function $f:U\setminus A\to\er$ has a unique extension on $U$. Is it then necessarily true that $A$ is intervally thin?
\end{problemTT}

Arguably our main result is that the answer to this question is negative. Example~\ref{HDS} and Remark~\ref{rem} provide a closed set $K$ which is not intervally thin, but which enjoys the ``unique extension property for convex functions'' (UEP) from Problem TT. We took the liberty of calling this set $K$ ``the Holey Devil's Staircase'' since it is the graph of the classical Cantor function (the Devil's Staircase) minus all the horizontal open line segments contained in the graph (in other words, it is the graph of the restriction of the Cantor function to the Cantor set).

One can readily verify that the Holey Devil's Staircase is not intervally thin. It is enough to consider the last intersection of the graph of the Cantor function with any line segment with endpoints in $(-\infty,0)\times\left(0,\frac12\right)$ and $(1,\infty)\times\left(\frac12,1\right)$; clearly, this intersection is an element of $K$.

To prove that $K$ has the UEP, is considerably more difficult and our effort in this direction has inspired a large part of this article.

The main result of \cite{Taborite} is essentially the following theorem. Note that since we restrict our attention to convex functions (as opposed to $\omega$-semiconvex functions studied in \cite{Taborite}), we change the formulation of the theorem accordingly:

\begin{thmTT}
Let $U$ be an open subset of $\er^d$ and let $A$ be a closed intervally thin subset of $U$. Let $f:U\setminus A\to \er$ be a locally convex function. Then $f$ has a unique locally convex extension on $U$.
\end{thmTT}

\noindent The proof of this theorem consists of two principal steps:
\begin{enumerate}[(1)]
\item First, one proves that there is a unique continuous extension; this is the more difficult part.
\item Once one has the continuous extension, it is then easy to prove that it is convex.
\end{enumerate}

Our aim is to apply this scheme to our set $K$. It turns out that in this case the easier step is (1); we only need a simple generalization of the corresponding theorem from \cite{Taborite}---which we have in Lemma \ref{ext}.

Performing step (2) for $K$ is the crucial part and it motivates the introduction of $c$-removable sets with the consequent natural question: Which sets are $c$-removable?

\begin{definition}
We say that a closed set $A\subset \er^d$ is \emph{$c$-removable} if the following is true: Every real function $f$ on $\er^d$ is convex whenever it is continuous on $\er^d$ and locally convex on $\er^d\setminus A$.
\end{definition}

A consequence of Theorem TT is that all closed intervally thin sets are $c$-removable, but this fact does not help us. In $\er^2$ we were able to find a sufficient condition more general than interval thinness which covers also the case of our set~$K$:

\begin{prop1}\label{prop1}
Let $K\subset\er^2$ be compact and intervally thin in two different directions.
Assume that for a dense set of line segments $L\subset \er^2$ the cardinality of
$K \cap L$ is at most countable. Then K is $c$-removable.
\end{prop1}
Here interval thinness of $K$ in a direction means that to any given line segment in that direction we can find arbitrarily close line segments contained in the complement of the set $K$. It is not difficult to see that for any closed set $K$ intervally thin in a direction $v$, any continuous function which is locally convex outside $K$ is necessarily convex on all lines parallel to $v$. Hence, the assumption of interval thinness of $K$ in two directions ensures that our function is convex in those two directions (i.e. is essentially separately convex) which we can use further in the proof---the key Lemma~\ref{sverak} tells us that a separately convex function cannot ``have a concave angle'' on any line.

The condition from the proposition may seem rather artificial, but it emerges quite naturally from our method of the proof. What is more, it is easily seen to be more general than interval thinness and is fulfilled by $K$. (Hence, the Holey Devil's Staircase is $c$-removable.) However, we were not able to generalize this condition to higher dimensions; instead, we used the geometric measure theory to obtain the following theorem which in $\er^2$ is strictly weaker than Proposition~1.

\begin{thm1}\label{theorem1}
Let $M\subset\er^d$ be a compact set which is intervally thin in $d$ linearly independent directions $n_1,\dots,n_d$.
Suppose that $M$ has $\sigma$-finite $(d-1)$-dimensional Hausdorff measure. Then $M$ is $c$-removable.
\end{thm1}

This condition does not include interval thinness because there are intervally thin sets of positive $d$-dimensional measure in $\er^d$. For instance, in $\er^2$ one can construct such a set by taking the full unit square and digging in it countably many straight tunnels in such a way that the rest is intervally thin but still of positive measure.

Among other signs, also from this fact it seems rather obvious that this theorem is far from being a characterization of $c$-removable sets. In fact, it is not even clear whether all $c$-removable sets in $\er^2$ are totally disconnected; from the considerations contained in the second part of Section~5 it seems plausible that the Koch curve might be an example of a non-trivial $c$-removable continuum in $\er^2$. (Of course, such an example has to be rather complicated as it is not difficult to prove that no smooth curve in $\er^2$ is $c$-removable.)

On the other hand, we have the following.

\begin{thm2}\label{theorem2}
Let $A, B\subseteq\er$ be closed sets of positive Lebesgue measure. Then $A\times B$ is not $c$-removable.
\end{thm2}
This theorem is interesting only for $A$, $B$ totally disconnected (otherwise $A\times B$ contains a non-degenerated line segment and the statement is trivial).
However, we do not know (and would like to know) whether e.g. the Cantor dust ($C\times C$ where $C$ is the Cantor set) is $c$-removable. As a matter of fact, possibly the most interesting of related open problems is:

\begin{problemLebNull}
Is there a closed totally disconnected Lebesgue null set in $\er^2$ which is not $c$-removable?
\end{problemLebNull}

It is worth pointing out that Theorem~2 is related to the recent work \cite{P} where a totally disconnected compact set which is not $c$-removable is constructed. The construction is rather complicated, but the witnessing function has a compact support, making the example stronger. However, even Theorem~2 is enough to achieve the main goal of \cite{P}, which is to disprove a theorem by L.~Pasqualini from 1938 \cite[Theorem~51]{Pasqualini} stating that any totally disconnected compact set in $\er^2$ is $c$-removable. It was the connection to this old article what convinced us that, of the two steps involved in the proof of Theorem~TT,  the crucial one is actually the second.

\section{Notation and basic facts}
All spaces shall be equipped with the Euclidean metric. We denote by $B(x,\eps)$ the open ball (with respect to the Euclidean metric) with the centre $x$ and radius $\eps$. Since confusion is unlikely, the symbol $(x,y)$ denotes an open interval in $\er$ as well as the point in $\er^2$ with coordinates $x$ and $y$. Similarly the symbol $[x,y]$ may denote a closed interval (when $x,y\in \er$) as well as the line segment with endpoints $x$ and $y$ (when $x,y\in \er^d$, $d>1$). By $\H^{k}$ we denote the $k$-dimensional Hausdorff measure.
For $M\subset\er^d$ and $\alpha$ a countable ordinal we denote the $\alpha$-th Cantor-Bendixson derivative of $M$ by $M^{(\alpha)}$.
The unit sphere in $\er^d$ is denoted by $S^{d-1}$.
For $v\in\er^d$ we denote the orthogonal complement of $v$ by $v^\bot$.
The symbol $\lin M$ denotes the linear span of $M\subset\er^d$.
For a fixed $d\in\en$ denote the standard basis of $\er^d$ by $\{e_1,\dots,e_d\}$.

Let $U\subset\er^d$ be open and $f:U\to\er$ be a function. We say that $f$ is \emph{locally convex} on $U$ if for some open convex $V\subset U$ the function $f|_V$ is convex. It is easy to see that a locally convex function is convex on any convex set contained in its domain.

The set $A\subset\er^d$ is called \emph{$k$-rectifiable} if there exist countably many Lipschitz mappings $f_i:\er^{k} \to \er^d$
such that
$$
\H^{k}\bigg(A\setminus \bigcup_{i=0}^\infty f_i\left(\er^{k}\right)\!\bigg)=0.
$$
Since we will work only with the case $k=d-1$, we will call $(d-1)$-rectifiable sets just \emph{rectifiable}.

Let $G(d,k)$ be the Grassmannian of $k$-dimensional linear subspaces of $\er^d$ equipped with the unique invariant probability measure $\nu_k^d$.
Besides the Hausdorff measure we will also use the $k$-dimensional Favard measure (integralgeometric measure) $\I^k$ on $\er^d$ which is for a Borel set $M$ defined as
$$
\I^k(M)=\frac{1}{\beta(d,k)}\int_{G(d,k)} \int _{V} \H^0(M\cap p_V^{-1}(y))\;d\H^{k}(y)\;d\nu_k^d(V),
$$
where $p_V$ is the orthogonal projection to $V$ and the number $\beta(d,k)$ is a non-zero constant depending only on $d$ and $k$ whose precise value is not important for us.

We will also need the following properties of the Favard measures.
Let $M\subset\er^d$ be a Borel set such that $\H^{d-1}(M)<\infty.$
Then $M$ can be expressed as a union of a rectifiable set $R$ and a set $P$ satisfying $\I^{d-1}(P)=0$ (c.f. \cite[3.3.13]{F}).
Moreover, each rectifiable set $R\subset\er^d$ satisfies $\I^{d-1}(R)=\H^{d-1}(R)$ (c.f. \cite[3.2.26]{F}).

\section{Separately convex functions}
The following lemma is a variant of an unpublished observation by V. \v Sver\'ak (see \cite{T}).
 For the convenience of the reader we provide a proof as we were not able to find one in the literature.
\begin{lemma}\label{sverak}
Let $f:\er^2\to\er$ be a separately convex function. Define $g:\er\to\er$ by $g(t)=f(t,t).$
Then
$$
 \liminf_{t\to 0+} \frac{g(x+t)+g(x-t)-2g(x)}{t}\geq 0
$$
for every $x.$
\end{lemma}

\begin{proof}
Without any loss of generality we can suppose that $x=(0,0)$, $f(z)=0$ and that there is a sequence $t_n\searrow0$ such that for each $n\in\en$,
\begin{equation}\label{direct}
 \frac{g(t_n)+g(-t_n)}{t_n}\leq -1.
\end{equation}
For $t>0$ put
$$
\sigma (t):=f(t,-t)+f(-t,t)
\quad\text{and}\quad
\rho (t):=f(t,t)+f(-t,-t).
$$
Note that, since $f$ is separately convex for every $t$,
\begin{equation}\label{zero}
\sigma(t)+\rho(t)\geq 0.
\end{equation}
Now we shall prove the following claim.

\emph{If for some $t_0$ and some $p$,
$$
\sigma(t_n)\geq pt_n \quad\text{for every $n$,}
$$
then
$$
\sigma(t_n)\geq (p+2)t_n \quad\text{for every $n$}.
$$
}
This is enough to prove the lemma since \eqref{zero} and \eqref{direct}, together with the above claim, imply $\sigma(t_n)\geq L$ for every $n$ and every $L\in\er$, which is not possible.

To prove the claim, first observe that due to \eqref{direct} we know that for each $n$,
$$
\rho(t_n)\leq -t_n.
$$
This implies for each $n$,
$$
\frac{\sigma(t_n)-\rho(t_n)}{2t_n}\geq \frac{p+1}{2}.
$$
By separate convexity of $f$, we get that
$$
f(t_k,t_n)+f(-t_k,-t_n)\geq \sigma(t_k)+ (t_n-t_k)\cdot\frac{\sigma(t_k)-\rho(t_k)}{2t_k}
$$
provided $k>n.$
If we consider $k\to\infty$, we obtain
$$
f(0,t_n)+f(0,-t_n)\geq t_n\frac{p+1}{2}.
$$
Using the separate convexity of $f$ one more time together with \eqref{direct} we get
\begin{equation*}
\sigma(t_n)\geq f(0,t_n)+f(0,-t_n)+(f(0,t_n)+f(0,-t_n)-\rho(t_n))\geq (p+2)t_n.
\qedhere
\end{equation*}
\end{proof}

\begin{definition}
The set $A\subset \er^d$ is called \emph{intervally thin in direction $v\in S^{d-1}$} if for any $x,y\in \er^d$ with $x-y$ parallel to $v$, and any $\eps >0$, there exist $x'\in B(x,\eps)$ and $y'\in B(y,\eps)$ such that $[x',y']\cap A=\emptyset$.
\end{definition}

\begin{proof}[Proof of Theorem~1]
First we claim that for every non-convex function $f$, there is a line $L$ such that $L\cap M$ is countable and $f|_L$ is non-convex. Moreover, $L$ can be found such that $v\in\lin\{n_i: i\in J\}$
for no $J\subsetneq\{1,\dots,d\}$, where $v$ is the direction of $L.$

To prove the claim first express $M$ as a countable union of sets $M_n$ satisfying $\H^{d-1}(M_n)<\infty.$
By \cite[3.3.13]{F} we can express each $M_n$ in the form $P_n\cup R_n$ with $\I^{d-1}(P_n)=0$ and
$R_n$ rectifiable.

Fix $n\in\en.$ Using \cite[3.2.26]{F} we see that
$$
\I^{d-1}(R_n)=\H^{d-1}(R_n)<\infty.
$$
This means, by the definition of the Favard measure, that for almost every $H\in G(d,d-1)$, almost every line perpendicular to $H$ intersects $R_n$ in at most finitely many points.
In particular, this means that almost every line intersects $M$ in finitely many points.

So, putting $P=\bigcup P_n$ and $R=\bigcup R_n$, we have that $\I^{d-1}(P)=0$, and also that almost every line intersects $R$ in at most countably many points.
Hence, almost every line intersects $M=P\cup R$ in at most countably many points.

Since $f$ is non-convex, the set $\A$ consisting of all lines $L$ such that $f|_L$ is non-convex has a positive measure.
To finish the proof of the claim we simply pick a line from $\A$ such that $L\cap M$ is at most countable.

Fix $z\in L$.
Due to the last part of the claim, we can suppose (possibly by composing $f$ with a suitable affine mapping) that $v_i=e_i$,
$v=\frac{1}{\sqrt d}(1,1,\dots,1)$  and such that $z=(0,0,\dots,0)$.

Now, we will prove the statement of the theorem using induction on $d$.
Suppose that $d=2$ and that $f$ is locally convex on $M^c$.
Put $K:=M\cap L$; then $K$ is a countable compact.
By Lemma~\ref{sverak} we know that if $f|_L$ is convex on $N^c$ for some $N$ then it is convex on a neighbourhood of any isolated point $N$.
This means that $f_L$ is convex on $K ^{(\alpha)}$ for every countable ordinal $\alpha.$
But, since $K$ is a countable compact, there is a countable ordinal $\beta$ such that $K^{(\beta)}=\emptyset$ which is a contradiction with the assumption of $f$ being non-convex.
This finishes the proof for $d=2.$

Suppose that the lemma is true for every $d$ up to $k-1\geq 2$, we will prove that it is true for $d=k$ as well.
Put $\nu:=v-v\cdot e_1$, $A(p):=pe_1+\lin\{e_2,\dots,e_d\}$ and $L(p):=pe_1+\lin\{\nu\}$ for $p\in\er.$
Then it is easy to verify that one of the following two statements is true:
\begin{enumerate}
\item[(a)] $f|_{L(p)}$ is non-convex for every $p$ from some interval $(a,b)$,
\item[(b)] $f|_{L(p)}$ is convex for every $p$.
\end{enumerate}
If $(a)$ is valid than by \cite[Theorem~7.7]{Mat} we know that $M\cap A(q)$ is of $\sigma$-finite $(d-2)$-dimensional Hausdorff measure for some $q\in(a,b)$.
But this is not possible by applying induction procedure to the function $f|_{A(q)}.$

One the other hand, $(b)$ is not possible either.
Indeed, we can apply Lemma~\ref{sverak} to $f|_{\lin\{v,e_1\}}$ the same way as in the proof of the case $d=2$.
\end{proof}

Note that the proof for $d=2$ directly gives us Proposition~\ref{prop1}.

\section{Extensions of locally convex functions}

\begin{definition}
We say that a set $A\subset\er^{d}$ is \emph{totally disconnected in a direction $v\in S^{d-1}$} if the set $A\cap l$ is totally disconnected for every line $l$ parallel to $v.$

\end{definition}

The following lemma is a refinement of \cite[Theorem 3.1]{Taborite} (note that the non-trivial part of the theorem is the existence of a unique continuous extension).

\begin{lemma}\label{ext}
Suppose that $A\subset\er^{d}$ is closed and both totally disconnected and intervally thin in some direction $v\in S^{d-1}$.
Let $U\subset\er^{d}$ be open.
Then every function locally convex on $A^{c}\cap U$ admits a continuous extension to $U.$
\end{lemma}

\begin{proof}
Let $f:U\to\er$ be locally convex on $A^{c}\cap U$.
Choose $x\in U$ and $\eps>0$ we need prove that there is a $\delta>0$ such that if $|x-a|,|x-b|\leq\delta$ and $a,b\in U\cap A^{c}$ then $|f(a)-f(b)|<\eps.$
Without any loss of generality we can suppose that $x=0$ and that $v$ is parallel to one of the coordinate axis.

For $u\in\er^{d}$ and $r>0$ put $l_u:=u+\lin\{v\}$ and $C(u,r):=u+[-r,r]^d$.
Since $A$ is totally disconnected in the direction $v$ we can find $\alpha>0$ and $\frac{\alpha}{2}>\gamma>0$ such that for $y:=\alpha v$ we have $C(y,\gamma),C(-y,\gamma)\subset A^{c}\cap U.$
Since $f$ is locally convex on $A^{c}$ and therefore locally Lipschitz on $A^{c}$, there is $K>0$ such that $f$ is K-Lipschitz on both $C(y,\gamma)$ and $C(-y,\gamma).$
Using the fact that $A$ is totally disconnected in the direction $v$ again, we can find $\min\left(\frac{\eps}{25K},\alpha-2\gamma\right)>\lambda>0$ such that for $z:=\lambda v$ we have $z\in A^{c}$.
Since $f$ is continuous on $A^{c}$, there is $\lambda>\delta>0$ such that for every $u\in C(z,\delta)$ we have $|f(z)-f(u)|\leq\frac{\eps}{4}$.

To obtain a contradiction, suppose that there are $a,b\in C(x,\delta)\cap A^{c}$ such that $|f(a)-f(b)|\geq\eps.$
Let $x_a$ and $x_b$ be the unique point in $(z+v^{\bot})\cap l_a$ and $(z+v^{\bot})\cap l_b,$ respectively.
Then $x_a,x_b\in C(z,\delta)$ and so $|f(z)-f(x_a)|\leq\frac{\eps}{4}$ and $|f(z)-f(x_b)|\leq\frac{\eps}{4}.$
Moreover, one of the inequalities
$$
f(a)-f(z)\geq\frac{\eps}{2},\;\;f(a)-f(z)\leq-\frac{\eps}{2},\;\; f(b)-f(z)\geq\frac{\eps}{2},\;\;f(b)-f(z)\leq-\frac{\eps}{2}
$$
must hold.
Therefore, one of the inequalities
\begin{equation}\label{nerovnice}
f(a)-f(x_a)\geq\frac{\eps}{4},\;\;f(a)-f(x_a)\leq-\frac{\eps}{4},\;\;f(b)-f(x_b)\geq\frac{\eps}{4},\;\;f(b)-f(x_b)\leq-\frac{\eps}{4}
\end{equation}
must hold as well.

Now, consider for instance the inequality $f(a)-f(x_a)\geq\frac{\eps}{4}$.
Since $A$ is intervally thin in the direction $v$ there are three colinear points $s_y\in C(-y,\gamma)$, $s_a\in C(x,\delta)$ and $s_x\in C(z,\delta)$ such that $[s_y,s_x]\subset U\setminus A$ and such that
$$
|f(s_a)-f(a)|, |f(s_x)-f(x_a)|\leq\frac{\eps}{16}.
$$
Then we have
\begin{equation}\label{eska}
|s_a-s_x|\leq |x-y|+2\delta=\lambda+2\delta\leq3\lambda\leq\frac{3\eps}{25K}.
\end{equation}
Moreover,
\begin{equation}\label{efka}
f(s_a)-f(s_x)\geq f(a)-f(x_a)-|f(s_x)-f(x_a)|-|f(s_a)-f(a)|\geq \frac{\eps}{4}-\frac{\eps}{16}-\frac{\eps}{16}=\frac{\eps}{8}.
\end{equation}
Using (\ref{eska}) and(\ref{efka}) we obtain
\begin{equation}\label{sklon}
\frac{f(s_a)-f(s_x)}{|s_a-s_x|}\geq \frac{\eps}{8}\cdot\frac{25K}{3\eps}>K.
\end{equation}
Choose an arbitrary $w\in([s_y,s_x]\setminus\{s_y\})\cap C(-y,\gamma)$.
From the convexity of $f$ on $[s_y,s_a]$, the fact that $[s_y,w]\subset[s_y,s_x]\subset[s_y,s_a]$ and \eqref{sklon} we obtain
$$
K<\frac{f(s_a)-f(s_x)}{|s_a-s_x|}\leq \frac{f(s_y)-f(w)}{|s_y-w|}\leq K,
$$
which is not possible. The remaining cases in \eqref{nerovnice} can be proved following the same lines.
\end{proof}

\begin{example}\label{HDS}
There is a compact set $K\subset\er^2$ which is not intervally thin and such that for every $f:K^{c}\to\er$ locally convex on $K^c$ there is a convex extension $F:\er^2\to\er.$
\end{example}
\begin{proof}
Let $h : [0, 1]\to [0, 1]$ be the classical Cantor function (the
Devil's Staircase) and let $C \subset[0,1]$ be the Cantor ternary set.
Now define the set $K$ as the graph of $h$ restricted to $C$.

First note that $\H^1(K)<\H^1(\graph h)<\infty$ and that $K$ is intervally thin in directions $(1,0)$ and $(0,1).$
Therefore using Theorem~\ref{theorem1} and Lemma~\ref{ext} we obtain that $K$ has the desired extension property and so it remains to prove that $K$ is not intervally thin.
Define $H:\er\to\er$ by $H=h$ on $[0,1]$, $H=0$ on $(-\infty,0)$ and $H=1$ on $(1,\infty).$
Then $\er^2\setminus\graph H$ has two components, say $C^+$ and $C^-$, and therefore for any $x^\pm\in C^\pm$ the line segment $[x^+,x^-]$ intersects $\graph H.$
Now take $x^+\in B((-\frac{1}{3},\frac{1}{3}),\frac{1}{4})$ and $x^-\in B((\frac{4}{3},\frac{2}{3}),\frac{1}{4})$ and set
$$
x:=\sup\{a\in\er: \text{there exists $b\in\er$ such that $(a,b)\in[x^+,x^-]\cap \graph H$}\}.
$$
Then $(x,H(x))\in K.$
\end{proof}

\begin{remark}\label{rem}
Note that the set $K$ from Example~\ref{HDS} also provides an answer to Problem~TT (see the introduction) in the negative.
This can be seen from the fact that the argument used in the proof of Example~\ref{HDS} can be easily localized using the  self affinity of $K$.
\end{remark}

\section{Two examples}
In the first part of this section we prove Theorem~2 which provides us with a very natural class of examples of sets which are not $c$-removable. We shall need the following definitions.

Let $\beta,\eps>0$. Then we define the function $g_{\beta,\eps}:\er^2\to\er$ by
$$
g_{\beta,\eps}(x,y):=
\begin{cases}
\beta y^2-2\eps x-\eps^2,\ &(x,y)\in(-\infty,-\eps]\times\er,\\
\beta y^2+x^2,\ &(x,y)\in(-\eps,\eps)\times\er,\\
\beta y^2+2\eps x-\eps^2,\ &(x,y)\in[\eps,\infty)\times\er.
\end{cases}
$$
For $w\in\er$, set
$$g_{\beta,\eps}^{w}(x,y):=g_{\beta,\eps}(x-w,y) \quad\text{and}\quad h_{\beta,\eps}^{w}(x,y):=g_{\beta,\eps}^{w}(y,x).$$
One can readily verify that all the functions just defined are convex and $C^1$.

The Hessian matrix of $g_{\beta,\eps}^{w}$ in $(-\infty,-\eps+w)\times\er$, $(-\eps+w,\eps+w)\times\er$ and $(\eps+w,\infty)\times\er$, respectively, is
\begin{equation}\label{hessvert}
\left(\begin{array}{cc}
	0 & 0\\
	0 & 2\beta
\end{array}\right),
\quad
\left(\begin{array}{cc}
	2 & 0\\
	0 & 2\beta
\end{array}\right)
\quad\text{and}\quad
\left(\begin{array}{cc}
	0 & 0\\
	0 & 2\beta
\end{array}\right),
\end{equation}
and similarly the Hessian matrix of $h_{\beta,\eps}^{w}$ in $\er\times(-\infty,-\eps+w)$, $\er\times(-\eps+w,\eps+w)$ and $\er\times(\eps+w,\infty)$, respectively, is

\begin{equation}\label{hesshor}
\left(\begin{array}{cc}
	2\beta & 0\\
	0 & 0
\end{array}\right),
\quad
\left(\begin{array}{cc}
	2\beta & 0\\
	0 & 2
\end{array}\right)
\quad\text{and}\quad
\left(\begin{array}{cc}
	2\beta & 0\\
	0 & 0
\end{array}\right).
\end{equation}

Further, define the functions $f_i:\er^2\to \er$, $i\in\{1,2,3,4\}$, as follows:
$$
f_1(x,y):=
\begin{cases}
\frac{1}{12}x^2+4(y-1)^2,\ &(x,y)\in\er\times (1,\infty),\\
\frac{1}{12}x^2,\ &(x,y)\in\er\times (-\infty, 1],
\end{cases}
$$
and $f_2(x,y):=f_1(x,-y)$, $f_3(x,y):=f_1(y,x)$ and $f_4(x,y):=f_3(-x,y)$.

Again, it is easy to check that the functions $f_i$, $i=1,2,3,4$ are $C^1$ and that the Hessian matrix of (e.g.) $f_1$ in $\er \times (1,\infty) $ and $\er \times(-\infty,1)$, respectively, is
\begin{equation}\label{hessvenku}
\left(\begin{array}{cc}
	\frac{1}{6} & 0\\
	0 					& 8
\end{array}\right)
\quad\text{and}\quad
\left(\begin{array}{cc}
	\frac{1}{6} & 0\\
	0 					& 0
\end{array}\right).
\end{equation}
It is also useful to note that the Hessian Matrix of $\varphi:(x,y)\mapsto -xy$ is
\begin{equation}\label{hessnekonv}
\left(\begin{array}{rr}
	 0 & -1\\
	-1 &  0
\end{array}\right).
\end{equation}

\begin{lemma}\label{LCRem}
Let there, for each $i\in\en$, be given $\eps_i>0$, $\beta_i \in \left(0, \frac{1}{80}\right)$ and $w_i\in(-1,1)$ such that $\sum\beta_i=\frac{1}{4}$, $\sum\eps_i<\frac{1}{24}$ and $(w_i-\eps_i,w_i+\eps_i)$ are pairwise disjoint intervals contained in $[-1,1]$.
Denote $g_i:=g_{\beta_i,\eps_i}^{w_i}$ and $h_i:=h_{\beta_i,\eps_i}^{w_i}$, define the function $f:\er^2\to\er$ as
$$
f(x,y):=\varphi(x,y)+\sum_{i=1}^4 f_i(x,y)+ \sum_{i=1}^\infty\left(g_i(x,y) +h_i(x,y)\right),
$$
and define the set $K\subset \er$ as
$$
K:=[-1,1]\setminus \bigcup_{i=1}^\infty (w_i-\eps_i,w_i+\eps_i).
$$
Then $f$ has the following properties:
\begin{enumerate}
\item It is non-convex, since
$$\frac{f(-1,-1)+f(1,1)}{2}<f(0,0);$$
\item it is locally convex on $\er^2 \setminus K^2$.
\end{enumerate}
\end{lemma}

\begin{proof}
First, we need to check that $f$ is a well-defined function. To that end, it is sufficient to note that for any $R>1$ and any $i\in\en$ the maximum of $g_i$ on $[-R,R]^2$ is attained at the point $(R,R)$ and is less than $\beta_i R^2+2\eps_i 2R$. Obviously, the same is true for $h_i$, and so the infinite series in the definition of $f$ converges locally uniformly on $\er^2$.

Now, the Hessian matrix of $f$ at a point $(x,y) \in \er^2\setminus K^2$ need not exist; however, all the summands in the definition of $f$, excluding $\varphi$, are convex functions. Since the sum of convex functions is convex, it is enough to prove that at any point $(x,y)$ outside $K^2$ we can find finitely many $f_i$'s, $g_i$'s and $h_i$'s such that their sum together with $-xy$ has a positively definite Hessian matrix at $(x,y)$.

This is easy to do if $(x,y)\notin[-1,1]^2$; in this case we only need to use one of the $f_i$'s. For example, if $(x,y)\in \er\times (1,\infty)$, then we see from \eqref{hessvenku} and \eqref{hessnekonv} that $\varphi+f_1$ is convex at $(x,y)$.

The ``worst case'' is that $(x,y)\in [-1,1]^2$ lies in a single vertical stripe of the form $(w_k-\eps_k,w_k+\eps_k)\times\er$ (fix the $k$) and not in any of the horizontal stripes of the form $\er\times (w_i-\eps_i, w_i+\eps_i)$.
In this case we find a finite set $F\subset\en$ such that $(x,y)$ does not lie on the boundary of any stripe of the form $\er\times (w_i-\eps_i, w_i+\eps_i)$ with $i\in F$ (which can happen at most twice) and such that
$$
\alpha_F:=\sum_{i\in F} \beta_i > \frac{9}{40}\,.
$$
Let us now consider the function
$$
f_F:=\varphi+\sum_{i\in F}\left(g_i+h_i \right) +
\begin{cases}
g_k, \ & \text{if } k\notin F; \\
0,   \ & \text{if } k\in F.
\end{cases}
$$
Since we know that $(x,y)$ does not lie on the boundary of any of the stripes (horizontal or vertical) involved in the definition of $f_F$, the Hessian matrix of $f_F$ exists on a convex open neighbourhood $U$ of $(x,y)$ and it is easy to see from \eqref{hessvert} and \eqref{hesshor} that its determinant satisfies
$$
\det\left(H_{f_F}(a,b)\right)\geq 4\left(\alpha_F+\alpha_F^2\right)-1>4\left(\frac{9}{40}+\frac{81}{1600}\right)-1>0, \quad (a,b)\in U.
$$
As $\frac{\partial^2 f_F}{\partial x^2}>0$, we obtain that $f_F$ is convex on $U$.
On the other hand, $f-f_F$ is convex and therefore $f$ is convex on $U$; one can check the other cases in a similar way concluding the proof of property $(b)$.

It remains to verify property $(a).$
Denote
$$
\delta_i:=\frac{g_i(-1,-1)+g_i(1,1)}{2}-g_i(0,0) .
$$
An easy computation shows that
$$
\delta_i=\beta_i+2\eps_i(1-|w_i|)
$$
and clearly the same is also true if we substitute all the occurrences of $g$ in the definition of $\delta_i$ by $h$. Hence,
$$
\sum_{i=1}^\infty \delta_i=\sum_{i=1}^\infty \beta_i+2\sum_{i=1}^\infty(\eps_i(1-|w_i|))<\frac{1}{4}+2\sum_{i=1}^\infty\eps_i<\frac{1}{3}\,.
$$
We also have
$$
\sum_{i=1}^4 \left(\frac{f_i(-1,-1)+f_i(1,1)}{2}-f_i(0,0)\right)=4\cdot\frac{1}{12}=\frac{1}{3}\,.
$$
The last two facts clearly imply $(a)$.
\end{proof}

\begin{proof}[Proof of Theorem~2]
If $A\times B$ contains a line segment then it is not $c$-removable as follows from \cite[Example 2.1]{Taborite}. From now on, assume that the sets $A$ and $B$ are totally disconnected.

By the Lebesgue density theorem we can find points $a\in \er$ and $b\in \er$ which are density points of $A$ and $B$ respectively. Without loss of generality assume that $a=b=0$; since $A$ and $B$ are both closed, $0\in A\cap B$. Then $0$ is clearly a point of density of $A\cap B\cap (-A)\cap(-B)$; consequently we can further assume that the sets $A$ and $B$ are symmetrical.  We shall prove that $(A\cap B)^2$ is not $c$-removable.

Take an $r>0$ such that $\frac{1}{2r}\lambda\left(A\cap B \cap (-r,r) \right) > \frac{23}{24}$ and such that $r$ (and therefore also $-r$) is in $A\cap B$. Without loss of generality we can assume that $r=1$. Now, for $i \in \en$ take $w_i\in(-1,1)$ and $\eps_i>0$ such that the intervals $(w_i-\eps_i,w_i+\eps_i)$ are pairwise disjoint and such that
$$[0,1]\setminus A\cap B= \bigcup_{i=1}^\infty(w_i-\eps_i, w_i+\eps_i).$$
The assumptions of Lemma \ref{LCRem} are now satisfied for any choice of positive numbers $\beta_i$, $i\in \en$, such that $\sum \beta_i = \frac{1}{4}$.
\end{proof}

The following two lemmas are concerned with the Koch curve and constitute a partial result regarding its $c$-removability. See also Problem~\ref{KochProblem}.

\begin{lemma}\label{limit}
$$
\lim_{k\to\infty}\left(3^k\prod\limits_{j=0}^{k-1}\frac{3^{j+1}+3}{3^{j+1}+1}
-2\sum\limits_{m=0}^{k-1}3^m\prod\limits_{j=0}^{m-1}\frac{3^{j+1}+3}{3^{j+1}+1}\right)=\infty\,.
$$
\end{lemma}

\begin{proof}
First, consider the following formula which follows easily by induction
$$
\prod\limits_{j=0}^{k-1}\frac{3^{j+1}+3}{3^{j+1}+1}=\frac{2\cdot 3^k}{3^k+1}\,.
$$
Now
$$
\begin{aligned}
2\left(\frac{3^{2k}}{3^k+1}-2\sum\limits_{m=0}^{k-1}\frac{3^{2m}}{3^m+1}\right)
=&2\left(3^k-\frac{3^{k}}{3^k+1}-2\sum\limits_{m=0}^{k-1}\left(3^m-\frac{3^{m}}{3^m+1}\right)\right)\\
=&2\left(3^k-2\sum\limits_{m=0}^{k-1}3^m-\frac{3^{k}}{3^k+1}+2\sum\limits_{m=0}^{k-1}\frac{3^{m}}{3^m+1}\right)\\
\geq&2\left(3^k-2\frac{3^k-1}{3-1}-1+2\sum\limits_{m=0}^{k-1}\frac{3^{m}}{3^m+1}\right)\\
=&2\left(1-1+2\sum\limits_{m=0}^{k-1}\frac{3^{m}}{3^m+1}\right)\geq 4\sum\limits_{m=0}^{k-1}\frac{1}{2}=2k.
\end{aligned}
$$
\end{proof}

\begin{lemma}
Suppose that $f$ is a continuous function locally convex on the complement of the Koch curve.
Then $f$ is convex on every line parallel to the $y$-axis.
\end{lemma}

\begin{proof}
First note that it is sufficient to prove the statement of the lemma for a dense set of lines parallel to the $y$-axis.
Due to the self similarity of the Koch curve it is then sufficient to prove that there is no continuous function on $[0,3]\times[0,\frac{6}{\sqrt{3}}]$
such that $f(\frac{3}{2},0)+f(\frac{3}{2},\frac{6}{\sqrt{3}})<2f(\frac{3}{2},\frac{3}{\sqrt{3}}).$
Here we consider the realisation of the Koch curve with endpoints $(0,0)$ and $(3,0)$.
For simplicity we will work in coordinates where the point $(\frac{3}{2},\frac{3}{\sqrt{3}})$ is translated to the origin.
For $i\in\en_0$ denote
$$
a_{i}=\left(0,-\frac{\sqrt{3}}{3^{i+1}}\right),\quad b_{i}=\left(\frac{1}{2\cdot 3^i},-\frac{\sqrt{3}}{3^{i+1}}\right),
\quad u_{i}=\left(\frac{1}{2\cdot 3^i},-\frac{\sqrt{3}}{3^{i+1}}\right),
$$
$$
z=\left(0,\frac{\sqrt{3}}{3}\right),\quad
s_{i}=\left(\frac{1}{2\cdot 3^i},\frac{\sqrt{3}}{3}\right)\quad \text{and}\quad
p_{i}=\left(\frac{1}{2\cdot 3^i},0\right).
$$
Modifying $f$ by adding an appropriate affine function and multiplying it by an appropriate constant we can suppose that
$f(z)\leq0$, $f(s_i)\leq0$, $f(a_0)=1$ and $f(p_i)\geq f(0,0)=1.$
Since $f$ is convex on $[s_i,b_i]$ for every $i$ we can write
\begin{equation}\label{one}
\begin{aligned}
f(b_i)\geq& f(u_{i+1})+\frac{|b_i-u_{i+1}|}{|s_i-u_{i+1}|}(f(u_{i+1})-f(s_{i}))\geq
\left(1+\frac{\frac{\sqrt{3}}{3^{i+1}}-\frac{\sqrt{3}}{3^{i+2}}}{\frac{\sqrt{3}}{3}+\frac{\sqrt{3}}{3^{i+2}}}\right)f(u_{i+1})\\
=&\left(1+\frac{2}{3^{i+1}+1}\right)f(u_{i+1})=\frac{3^{i+1}+3}{3^{i+1}+1}f(u_{i+1}).
\end{aligned}
\end{equation}
Moreover, since $f$ is convex on $[a_i,u_i]$ for every $i$ we can write
\begin{equation}\label{two}
\begin{aligned}
f(u_i)\geq& f(b_{i})+\frac{|b_i-u_{i}|}{|b_i-a_{i}|}(f(b_{i})-f(a_{i}))\\
\geq&f(b_i)+\left(\frac{\frac{1}{2\cdot 3^i}-\frac{1}{2\cdot 3^{i+1}}}{\frac{1}{2\cdot 3^i}}\right)(f(b_{i})-1)=3f(b_i)-2.
\end{aligned}
\end{equation}
Combining \eqref{one} and \eqref{two} we then obtain
$$
f(b_i)\geq \frac{3^{i+1}+3}{3^{i+1}+1}(3f(b_i)-2)\,,
$$
and iterating for every $i<k$,
$$
f(b_i)\geq 3^k\prod_{j=0}^{k-1}\frac{3^{j+1}+3}{3^{j+1}+1}f(b_{i+k})
-2\sum\limits_{m=0}^{k-1}3^m\prod\limits_{j=0}^{m-1}\frac{3^{j+1}+3}{3^{j+1}+1}\,.
$$
Since $f$ is convex on $[s_i,b_i]$ we have for every $i$,
$$
f(b_i)\geq f(p_{i+1})+\frac{|s_i-p_{i+1}|}{|b_i-p_{i+1}|}(f(p_{i+1})-f(s_{i}))\geq f(p_{i+1})\geq 1.
$$
Finally, for $i=0$  and any $k>0$ we obtain
$$
f(b_0)\geq 3^k\prod\limits_{j=0}^{k-1}\frac{3^{j+1}+3}{3^{j+1}+1}
-2\sum\limits_{m=0}^{k-1}3^m\prod\limits_{j=0}^{m-1}\frac{3^{j+1}+3}{3^{j+1}+1}
$$
which is not possible due to Lemma~\ref{limit}
\end{proof}

\section{Open problems}

\noindent The following general question is likely to be very difficult to answer, but naturally arises from the introducion of the notion of $c$-removability.

\begin{problem} \label{ProbUltimatni}
Is there any interesting characterization of $c$-removable sets?
\end{problem}

However, there are several other interesting problems whose solutions might contribute to our understanding of the matter.

\begin{problem}\label{NullProblem}
Is there a closed totally disconnected Lebesgue null set in $\er^2$ which is not $c$-removable?
\end{problem}

\begin{problem}
Is the Cantor dust $c$-removable?
\end{problem}

\begin{problem}\label{KochProblem}
Is there a non-trivial $c$-removable continuum in $\er^2$?
\end{problem}

Note that if one could prove that there is a dense set of lines intersecting the Koch curve in countably many points, the answer would be positive; this would follow from the proof of Theorem~\ref{theorem1}.


\section{Introduction}
The notion of scale-free networks has been around for a
while~\cite{dP76}. The introduction of the
\emph{preferential attachment} (PA) algorithm for generating random
scale-free graphs was a significant step in understanding the properties
of scale-free networks, and the physical processes that create them \cite{aB99}. PA
has spawned a good deal of subsequent algorithms and  analysis.
The purpose of this paper is to highlight the straightforward fact that not all scale-free networks
have the same properties, and that algorithms, like PA  (for example), do not capture the
full richness of scale-free networks, nor do they necessarily display
properties that may be termed typical of all scale-free networks.

To achieve our aim, we first briefly recall in this introduction the
principal processes that have been proposed to describe and generate
scale-free networks, and indicate deficiencies of these processes when employed as models of
typical scale-free networks. We then propose a \emph{maximum entropy
  process} that provides an unbiased sample of the set of all scale-free
networks. A \emph{maximum entropy process} provides a better
representation of expected properties of scale-free networks, both in
terms of richness and typicality. A \emph{maximum entropy process} also
provides a unbiased standard against which other processes that generate
scale-free networks can be compared. While a great deal of this maximum entropy process reduces to simple edge-switching, it is the application of this process to achieve maximum entropy realisations that is important.

In Sec.~\ref{sec:results} we make a careful comparison of PA with our
unbiased standard to illustrate how PA has a
significant bias in the structural properties of the networks it
generates. We demonstrate it has a \emph{hub-centric} bias.
In Sec.~\ref{sec:applications} we use a surrogate data approach to examine a
selection of real-world networks claimed to be scale-free, and analyze in
what sense these networks are typical of scale-free networks and how they
differ.

\subsection{Scale-free networks from processes}

In this section we briefly review the notion of scale-free networks, and
the principal models of the physical processes that generate scale-free
networks. These models can be broadly divided into \emph{growth} models
and \emph{configuration} models.

Scale-free networks are usually identified by the histogram of node
degrees having a power-law tail. Many naturally occurring networks{\footnote{For a sample of the data we use here, see {\tt http://vlado.fmf.uni-lj.si/pub/networks/data/} .}} have
been identified as having a scale-free property: citation and
collaboration networks~\cite{dP65,aW04}, airline
networks~\cite{rG04}, protein-protein interaction~\cite{hJ01}, metabolic
pathways~\cite{jD05}, the world-wide web and internet~\cite{aH13}.

To understand the formation of scale-free networks various models have
been proposed to mimic the physical or conceptual processes that build and
shape these networks. One of the first, proposed by Barab\'asi and Albert,
is \emph{preferential attachment}~\cite{aB99}, which is a
restatement of the process described by de Solla Price~\cite{dP76} as
a model of observed scale-free networks of citations~\cite{dP65}.

Preferential attachment is an \emph{unchanging}, \emph{additive} growth
process, where nodes of a fixed degree are added to the network, with
links preferentially attached to existing nodes depending on their degree;
usually proportional to the degree. Although this was suggested as a model
of the process that grows the world-wide web (seen as hyperlinks between
webpages), it was recognized that there is an additional \emph{aging}
process (AOL is replaced by Facebook, yahoo looses popularity to Google)
so that attachment preferences \emph{change} as the network growth
{proceeds}~\cite{rX04,pK01,dP02}.

Other processes can be introduced into a network growth model, such as
\emph{shuffling} parts of the network, and \emph{deletion} of links and
nodes~\cite{rX04, kJ13}, such that these actions do not change
the underlying scale-free property of the networks produced. Moreover,
scale-free networks can be produced by processes that
are not growth processes. The most commonly considered of these are the so-called {\em configuration models}.  Configuration models proceed by choosing nodes
to have prescribed degrees, then connecting these together to form a
scale-free network~\cite{mM95}; although care is needed to ensure
the networks obtained satisfy other expected properties, such as, being
simple (no self-links or multiple edges between nodes) and
connected~\cite{tB06, pC12, kJ13}. In addition to the problems with multplie edges and self-links \cite{kJ13}, configuration models  do not provide a well-founded sampling --- in the sense of achieving maximum entropy realisations \cite{gB09}.

\subsection{Not all scale-free networks are the same}

With so many different processes generating scale-free networks, the
question arises: do they generate the same type of networks? The simple answer is
no. Many differences in the scale-networks generated by different
processes have been noted; some particular differences are outlined in the
following.

Consider a preferential attachment process ${\rm PA}(m,A)${} that {{attaches}} nodes of
degree~$m$ with the probability of linking a new node to a node of
degree~$k$ is $k+A-1$, where $A\geq1$ is a constant called the
\emph{initial attractiveness}~\cite{sD00c}. The minimum degree~$m$
clearly effects the \emph{robustness} of the networks generated, and it
effects other properties too~\cite{kJ13}. It can be shown that for ${\rm PA}(m,A)${}
the power-law tail has an exponent $\gamma=2+A/m$. This type of
preferential attachment process restricts the power-law exponent to
$\gamma>2$. Many real-world networks exhibit $\gamma<2$. In such cases the graphs have been termed {\em dense}~\cite{cG11} and generating such graphs has been seen as problematic. This led Del Genio and colleagues to conclude eponymously  to the contrary.  However, we recently provided an algorithm that easily find dense graph --- the apparent paradox with the claim of \cite{cG11} is resolved when one observes that \cite{kJ13} searches over the space of all networks while \cite{cG11} starts with the constraint of viable graphs of a fixed degree sequence. Hence, both \cite{fF13} and \cite{kJ13} proposed that these network growth models can be modified to shape networks and allow such
power-laws.

When other network statistics are considered (such as diameter, node-degree
assortativity, clustering coefficients), then the differences between
scale-free networks found in the real-world and generated by different
models becomes more apparent.

For example, consider the node-degree assortativity of scale-free networks
generated by preferential attachment. These networks have low
assortativity, but many natural networks are found to have very high
assortativity~\cite{mS08}. Although these natural networks are not
small-world (although they are scale-free typically estimates of the exponent are much smaller than observed for most processes that generate scale-free networks:  $\gamma<2$), it can be shown that other processes can generate scale-free
networks with high node-degree
assortativity~\cite{skinny2}. Newman~\cite{mN03} observed systematic
biases in the assortativity of real-world scale-free networks:
technological networks tended to be dis-assortative, social networks
assortative. He also showed that networks generated by preferential
attachment have a definite bias compared to real-world networks; on the
other hand, \emph{altruistic attachment} does not~\cite{pL11}. Other issues
with the assortativity of preferential attachment networks ~\cite{rX04,xu3} and scale-free networks more generally \cite{sJ10} have been
observed previously.

It has been noted that growth models have systematic
biases (of network statistics) relative to configuration
models~\cite{dC01,sD00c}.  However, while configuration models can successfully alleviate that bias, it is not clear that they provide a sampling of the appropriate distribution. In particular, networks of growth models
usually only attain the scale-free property asymptotically, and small
networks display systematic biases~\cite{xu3,mS08};  configuration
models remove this particular growth-based bias, but it is not clear at what expense. For this reason we seek a maximum entropy realisation from the ensemble of viable graphs \cite{gB08}.

\subsection{Maximum entropy processes}

In this section we aim to define a process that generates scale-free
networks in as pure form as possible, that is, without any of these biases. This
process is a \emph{maximum entropy} process, which can serve as a standard
to which other processes can be compared to both recognize and understand
the nature of their particular biases.

As previously noted, \emph{scale-free} is usually taken to mean that the
node-degree histogram has a power-law tail. We need to make this
definition more precise. Let ${\mathcal{G}}$ be the set of connected simple graphs
with $N$ nodes, and let ${\mathbf{n}}(G)=(n_1,\dots,n_{N-1})$, $n_k\in\Z$, be the
degree histogram of~$G\in{\mathcal{G}}$, that is, $n_k$ is the number of nodes of
degree~$k$.

\begin{defn}\label{defn:tail}
  If ${\mathbf{p}}=(p_1,\dots,p_{N-1})$, where $p_k\in[0,1]$ and
  $\sum_{k=1}^{N-1}p_k=1$, then ${\mathbf{p}}$ has a \emph{power-law tail} if
  $p_k\propto(k-\alpha)^{-\gamma}$ for $k\geq\beta$, where
  $\alpha,\beta\in\Z$, $\gamma\in\R$, $\alpha,\beta\geq0$ and
  $\gamma>1$.
\end{defn}

In this definition $\alpha$ shifts the power-law
tail relative to~$k$,  $\beta$ is where
the power-law tail is deemed to begin, and
$\gamma$ is the power-law of the tail.
For the special case of preferential attachment ${\rm PA}(m,A)${}, the key parameters
$\alpha$, $\beta$ and $\gamma$ are implicitly defined by $m$ and~$A$.

\begin{defn}\label{defn:sfgraph}
  $G\in{\mathcal{G}}$ is a \emph{scale-free graph} if ${\mathbf{n}}(G)\simeq N{\mathbf{p}}$, where ${\mathbf{p}}$
  has a power-law tail.
\end{defn}

Although Defn.~\ref{defn:sfgraph} is apparently what most researchers appear to mean by \emph{scale-free}, what \emph{approximately equal} means
in~${\mathbf{n}}(G)\simeq N{\mathbf{p}}$ is not always clear. It is usually taken to mean that
the form $(k-\alpha)^{-\gamma}$ \emph{fits} ${\mathbf{n}}(G)$, in the sense of
least-squares curve-fitting of $\log{}n_k$ against $\log{}k$, or a
multinomial fit of ${\mathbf{n}}$ against~${\mathbf{p}}$.

Of interest are processes generating independent random scale-free
graphs. A stationary stochastic process~${\mathcal{P}}$ generating independent random
graphs in~${\mathcal{G}}$ is equivalent to selecting graphs according to a fixed
probability mass function ${\mathcal{P}}:{\mathcal{G}}\to[0,1]$, where
$\sum_{G\in{\mathcal{G}}}{\mathcal{P}}(G)=1$. Write $G\sim{\mathcal{P}}$ to denote that $G\in{\mathcal{G}}$ is
selected according to process~${\mathcal{P}}$. (Since the process is completely
defined by its probability mass function, we use the same symbol~${\mathcal{P}}$ to
denote the process and its probability mass function.)

\begin{defn}\label{defn:sfprocess}
  A process~${\mathcal{P}}$ \emph{generates} scale-free graphs (on average), if
  \begin{equation}
    \label{eq:sf}
    \E[\,{\mathbf{n}}(G) \mid G\sim{\mathcal{P}}\,] = N{\mathbf{p}},
  \end{equation}
  where ${\mathbf{p}}$ has a power-law tail.
\end{defn}

Such a process generates \emph{scale-free graphs}, according to
Defn.~\ref{defn:sfgraph}, with the notion of \emph{approximately equal},
${\mathbf{n}}(G)\simeq{}N{\mathbf{p}}$, implicitly defined by~${\mathcal{P}}$. For example, one can
compute the variance
\[\E\left[\,\sum_{k=1}^{N-1}(n_k-Np_k)^2 \Biggm| G\sim{\mathcal{P}}\,.\right],\]
or compute the probability of a histogram~${\mathbf{n}}(G)$ deviating from $N{\mathbf{p}}$ in
some norm by more than a threshold~$\eps$, that is,
$\Pr(\|{\mathbf{n}}(G)-N{\mathbf{p}}\|>\eps \mid G\sim{\mathcal{P}})$.

Preferential attachment ${\rm PA}(m,A)${}  is one example of such a process; a network is
grown by random preferential attachment until a graph with $N$ nodes is
obtained. This process implicitly defines the probability mass function
for $G\in{\mathcal{G}}$, which depends on~$m$ and on how the seed network is
chosen. Other growth models (with link ageing, initial attractiveness,
shuffling) will usually obtain different probability mass functions.

A difficulty with growth processes, including ${\rm PA}(m,A)${}, is that although they
satisfy Defn.~\ref{defn:sfprocess} of scale-free processes, the
power-law~$\gamma$ is implicitly defined. In some cases (although not for
PA) it is not clear that $\gamma$ is necessarily the same for all
graphs generated for fixed growth parameters. Depending on the initial
seed graph, it is possible to imagine growth processes which would
converge to different power-laws. In contrast, one of the advantages of
configuration models is the mean histogram~$N{\mathbf{p}}$ is prescribed: the
process configures nodes with pre-defined degrees into a
network~\cite{mM95}.

Even if two processes generate scale-free graphs by
Defn.~\ref{defn:sfprocess}, or have the same power-law~$\gamma$, or have
the same expected histogram~$N{\mathbf{p}}$, the expected value of other graph
statistics such as assortativity, clustering, diameter, may be quite
different --- as noted above. If we are to understand the nature of scale-free graphs, and the nature of
the processes that generate them, then it is useful to establish some
\emph{benchmark} or \emph{standard} against which processes are
compared. To some extent ${\rm PA}(m,A)${} has served that propose, principally
because it is a simple growth process. Unfortunately, it is not a
satisfactory standard for many reasons we already indicated: it has a
limited range of~$\gamma$ ({{\em typically} $\gamma$ is constrained to be between $2$ and $3$}; for $\alpha\gg m$ and $N\rightarrow\infty$ one obtains $\gamma\rightarrow 3$); ${\mathbf{p}}$ and $\gamma$ are implicitly defined; and, is subject to
well known biases of several basic network measures (such as assortativity).

We propose that the best standard against which to compare processes that
generate scale-free networks is a \emph{maximum entropy}
process~${\rm ME}(\p)${} with prescribed mean histogram~$N{\mathbf{p}}$. A maximum entropy
process samples all graphs with mean histogram~$N{\mathbf{p}}$ with the least bias
or emphasis of structures and features that are not common features of all
such graphs.

\begin{defn}\label{defn:maxentropy}
  The \emph{maximum entropy process}~${\rm ME}(\p)${} uniformly samples ${\mathbf{n}}$
  by the multinomial distribution of~${\mathbf{p}}$, and uniformly samples graphs
  with the histogram~${\mathbf{n}}$. That is, $\Pr(G \mid G\sim{\mathcal{P}}_{{\rm ME}({\mathbf{p}})}) =
  \Pr(G\mid{\mathbf{n}})\Pr({\mathbf{n}}\mid{\mathbf{p}})$, where
  \begin{equation}
    \label{eq:uniform}
    \Pr(G\mid{\mathbf{n}}) = \frac{1}{\left|\left\{G\in{\mathcal{G}}\mid{}{\mathbf{n}}(G)={\mathbf{n}}\right\}\right|}
  \end{equation}
  and
  \begin{equation}
    \label{eq:multinomial}
    \Pr( {\mathbf{n}} \mid {\mathbf{p}}) = N!\prod_{k=1}^{N-1}\frac{p_k^{n_k}}{n_k!}.
  \end{equation}
\end{defn}
The process~${\rm ME}(\p)${} generates connected graphs with maximum entropy from the specified degree sequence
partitions all graphs in~${\mathcal{G}}$ into equivalence classes by histogram~${\mathbf{n}}$,
and uniformly samples these~(\ref{eq:uniform}), hence, achieving maximum
entropy within these equivalence classes.  Moreover, the sampling of~${\mathbf{n}}$ uses a
multinomial distribution~(\ref{eq:multinomial}), which is the maximum
entropy process for selecting a histogram with frequencies~${\mathbf{p}}$. Hence, even
though the degrees of nodes within a graph are not independent, the
degrees of nodes of successive randomly generated graphs behave as though
they were independent. That is, maximum entropy is a consequence of selecting degree histogram by multinomial sampling and thence uniform sampling from within that degree histogram. This is equivalent to the work presented by Bianconi \cite{gB08, gB09} for selecting maximum entropy for a specified ensemble. For us, degree distribution and connectivity constraints dictate the ensemble of interest.

For details on how to efficiently implement~${ME({\mathbf{p}})}$, see
Algorithm~\ref{alg:MEp} in \ref{App:mainalg}. A \emph{maximum
  entropy} process~${ME({\mathbf{p}})}$ can be implemented by as a Monte Carlo
Markov Chain (MCMC)~\cite{kJ13}, however, a naive MCMC implementation
could converge slowly and not be effective for large
graphs~\cite{pD12, pD13, yC05}. The key ideas of
Algorithm~\ref{alg:MEp} are that the equivalence class of graphs with the
same degree histogram can be sampled uniformily by a process of \emph{edge
  switching}~\cite{aS11, rT81}, and a cannonical graph with a
given degree histogram can be constructed by the Havel-Hakimi
process~\cite{vH55, sH62}. Combining these key ingredients with a test
for connectivity~\cite{kJ13} provides a sufficient algorithm.

Following Bianconi \cite{gB08,gB09} we define the maximum entropy scale free network ensemble as follows.
\begin{defn}\label{defn:UniS}
  Let $SF(\gamma,m)$ denote the maximum entropy process $ME({\mathbf{p}})$, where
  ${\mathbf{p}}$ has a power-law tail~$\gamma$, with $\alpha=\beta=m$ and $p_k=0$
  for $k<m$.
\end{defn}
In the following $SF(\gamma,m)$ will serve as our standard against which
scale-free graphs and scale-free graph generators are compared.

\section{A hub of hubs is more than a rich club}
\label{sec:results}

In this section we use algorithm~\ref{alg:MEp} to explore what are the
typical properties of scale-free networks, and moreover, infer that PA
networks exhibit an atypical ``hub-centric'' structure. Power-law
distributions alone imply that the degree distribution has a long tail and
some nodes in the tail must have high degree, but the presence of ``hubs''
implies something more and are often said to ``hold the network
together''. What we observe in the following is something more than
this. There are three properties to which we refer with the term
``hub-centric'' (which we define later in an algorithmic sense): (i) the
distribution of hubs throughout the network, (ii) their interconnection,
and (iii) their connection with low-degree nodes. The presence or absence
of these three properties determine, to a very great extent, many of the
global properties of a scale-free network. Moreover, we observe that PA
networks exhibit properties which are atypical of the broader class.

However, before moving to the conclusion that the difference between these
two types of networks is due to the hub-centric properties of PA networks,
we first present a numerical study of the various main measures of network
geometry: assortativity, clustering coefficient and diameter. We apply
these measures to various different families of scale-free networks to
highlight the wide variation which is possible. We will also explore other
prominent features --- most importantly: (i) the robustness to targeted
worst-case attacks (i.e. ``attack vulnerability''
\cite{aR00}) on hubs; and (ii) a detailed motif
analysis.

For our comparison experiments we generate four types of PA$(m,A)$
networks, with $m=1,2,3,4$ and $A=0$. For given $N$ and~$m$ we generate a
sample from PA$(m,A)$ of~$M$ graphs~$G_i$, $i=1,\dots,M$. From a degree
histogram~${\mathbf{n}}(G)$ we use the following formula~\cite{aC09} to
estimate the power-law tail exponent~$\gamma$ of~$G$, as
\begin{equation} \label{eq:est_g} \widehat\gamma(G) = 1+N
    \left(\sum_{k=1}^{N-1} n_k\log\frac{k}{m-\frac12}\right)^{-1}.
 \end{equation}
Letting $\gamma=\frac1M\sum_{i=1}^M\widehat\gamma(G_i)$ we use
Algorithm~\ref{alg:MEp} and generate an sample of~$M$ graphs from
${\rm SF}(\gamma,m)${}. Results reported in the following are for $N=2000$, $m=1$, $2$, $3$,
$4$, and $M=40$. Computation with other values produced comparable
results. In the following subsections we examine each of the network
properties outlined above, starting with motif rank and then robustness.

\subsection{Comparison of motif rank }

The motif, defined as a small connected subgraph that recurs in a graph, is the basic building block, or functional unit, of complex networks \cite{aJ09}. In real-world networks (e.g. gene regulatory networks), the relative frequencies of different motifs often represent different functions of the network \cite{rM02,uA07}. Moreover, it sheds light on the structure of the whole network. We restrict our analysis here to the four-node motifs, and classify them into three groups (based on the number of edges in the constituent motifs, see  Fig.~\ref{motif}): $\{A,D\}$, $\{B,E\}$, $\{C,F\}$. Measuring the relative frequency of each of these groups allows us to determine the proportion of building blocks typical of the sparse, moderate, and dense networks. The result of the motif analysis is also shown in Fig.~\ref{motif}.

\begin{figure}
\begin{center}
\includegraphics[width=0.95\textwidth]{figs/fig1.png}
\end{center}
\caption{Each figure shows the frequencies of six types of four-node motifs, for a  corresponding networks with $\textbf{m}\in\{1,2,3,4\}$. On the left, are the results for preferential attachment, on the right, the uniform sampling algorithm.  For each network, ${\rm PA}(m,A)${} or ${\rm SF}(\gamma,m)${} (denoted in the figure as UniS) with specified $m$, we compute the frequency of each of the size possible arrangements of four-node path-connected subgraphs (labelled A to F) and then plot the absolute (logarithmic) frequency of occurrence of each graph. The order from most to least frequent  is also provided.}
\label{motif}
\end{figure}

From Fig.\ref{motif}, there is an apparent difference between the motif distribution of these two type of networks: as $m$ increases, the motif frequencies of PA networks become larger, and the ``dense'' motifs of type $C$ occur more frequently than the less dense type $E$ --- suggesting that the PA networks are denser than the corresponding ${\rm SF}(\gamma,m)${} networks. This would be a natural consequence of the PA networks becomes more hub-centric (more cross links between hubs) as $m$ increases, while the hubs of the ${\rm SF}(\gamma,m)${} network remain more evenly distributed (less cross-links between hubs, and correspondingly fewer ``dense'' motifs).

\subsection{Robustness and attack vulnerability}

Due to the long tail of the power-law distribution, scale-free networks
are often claimed to be simultaneously robust to the random loss of nodes
(i.e. ``error tolerance'') and fragile to targeted worst-case attacks
\cite{aR00}. The robustness is seen to be a
consequence of the extremely high number of (``unimportant'') low degree
nodes, the fragility is due to the extremely important role of the
hubs. { This property is also called the ``Achilles' heel'' of PA
networks, or the ``robust yet fragile'' feature}. Intuitively, the
inhomogeneous connectivity distribution of many networks caused by the
power-law distribution would possess these two properties. However, from
our analysis, ${\rm SF}(\gamma,m)${} networks generated by sampling uniformly from the
family of all scale free networks do not exhibit this second
property. Again, we may attribute this to the hub-centric nature of PA
networks.

We quantify the robustness to targeted attacks by specific removal of the most highly connected (or important) nodes until the network is disconnected. We then take the number of the nodes removed as a measure. For illustration, we use degree and betweenness centrality (BC) of a vertex as measures with which to target nodes for removal. Roughly, BC is defined as the number of geodesics (shortest paths)
going through a vertex:
\[{\rm BC}(v)= \sum _{ i\neq j,i\neq v,j\neq v }^{  }{ { g }_{ ivj }/{ g }_{ ij } },\] where $g _{ ij } $ is total number of shortest paths between node $i$ and $j$, and $g _{ ivj } $ is the number of those paths that pass through $v$. The results are plotted in  Fig.~\ref{robu}.

The case of targeted attack is trivial for $m=1$, in which situation the networks are highly likely to be trees. Hence, we restrict our analysis for the case where $m=2,3$. these results are shown in Fig.~\ref{robu}.
From Fig.~\ref{robu}, it is safe to conclude that the ${\rm SF}(\gamma,m)${} network is much more robust than the corresponding PA network when facing targeted attack. This is a consequence of the fact that the PA networks are more hub-centric, while the ${\rm SF}(\gamma,m)${} networks are not. The fragility of the PA networks is due to the placement of hubs within the network --- not the scale-free-ness of the network {\em per se.}

\begin{figure}[!htb]
\begin{center}
\includegraphics[width=0.75\textwidth]{figs/fig2.pdf}
\end{center}
\caption{Numerical estimates of histograms showing the number of nodes which must be removed to induce collapse of the network giant component. Nodes are selected for removal via either node degree or node betweenness centrality (here we illustrate results for removal via node degree --- results for betweenness centrality are visibly indistinguishable). Panels (from top to bottom) are for different minimum degrees $m=1,2,3$. Results for PA are plotted as a solid blue line and and results for ${\rm SF}(\gamma,m)${} are plotted as a red dashed line. Histograms are computed with logarithmically equally space bins --- the same binning is used for both network construction methods.
We find that, ${\rm SF}(\gamma,m)${} networks are, for $m>1$, an order of magnitude more robust to targeted attack than typical PA networks.}
\label{robu}
\end{figure}

\subsection{Numerical statistics of network}

In this section we present the difference between the PA networks and ${\rm SF}(\gamma,m)${} networks by computing the widely quoted numerical network properties:
\begin{itemize}
\setlength{\itemsep}{0pt}
\setlength{\parsep}{0pt}
\setlength{\parskip}{0pt}
\item \emph{diameter}: the maximal shortest path length.
\item \emph{global clustering coefficient}: the proportion of the number of closed triplets divided by the number of connected triples of vertices, which is a measure of degree to which nodes in a graph tend to cluster together.
\item \emph{local clustering coefficient}: the proportion of links between the vertices within its neighborhood divided by the number of links that could possibly exist between them, which quantifies how close its neighbors are to being a clique.
\item \emph{assortativity}: the Pearson correlation coefficient of degree between pairs of linked nodes, which measures the preference for a network's nodes to attach to others that are similar in degree.
\end{itemize}
The application of these four statistics is summarized in Fig.~\ref{box}.

\begin{figure}
\begin{center}
\includegraphics[width=0.95\textwidth]{figs/fig3.pdf}
\end{center}
\caption{Boxplot of the four network statistics described in the text { and Table \ref{tab:tab1}}. The boxplot depicts maximum and minimum ({extremum of the box), upper quantile ($75\%$) and lower quantile ($25\%$) (the darker box), and median (solid dark line within the inner box)} of the data. The { left-most box within} each pair is the summary of the ${\rm SF}(\gamma,m)${} network, while the right one represents PA networks.  Computation is depicted for $m=1,2,3,4$. Note that, for example, assortativity has an increasing negative bias for PA networks, while ${\rm SF}(\gamma,m)${} approaches $0$ as $m$ is increased.}
\label{box}
\end{figure}

\begin{sidewaystable}
  \centering
    \begin{tabular}{r|r|r|r|r|r}
    \toprule
          & diameter & global\_cc & local\_cc & assortativty & motif rank \\
    \midrule
    m=1 PA & 19.47 (1.71) & 0     & 0.67 (0.01) & -0.095 (0.017) & BECFDA \\
    m=1 ${\rm SF}(\gamma,m)${} & 15.55 (4.33) & 0.0037 (0.0033) & 0.73 (0.03) &  -0.18 (0.05) & FECBDA \\
    m=2 PA & 7.05 (0.22) & 0.0059 (0.0005) & 0.061 (0.011) & -0.14 (0.014) & FCEBDA \\
    m=2 ${\rm SF}(\gamma,m)${} & 9.45 (1.53) & 0.014 (0.004) & 0.093 (0.14) & -0.095 (0.052) & FECBDA \\
    m=3 PA & 5.28 (0.45) & 0.011 (0.0004) & 0.12 (0.014) & -0.19 (0.015) & FECBDA \\
    m=3 ${\rm SF}(\gamma,m)${} & 7.57 (0.54) & 0.0087 (0.0021) & 0.017 (0.014) & -0.038 (0.019) & FCEBDA \\
    m=4 PA & 4.93 (0.27) & 0.014 (0.0004) & 0.18 (0.017) & -0.24 (0.019) & FECBDA \\
    m=4 ${\rm SF}(\gamma,m)${} & 6.67 (0.47) & 0.0072 (0.0019) & 0.011 (0.011) & -0.021 (0.017) & FCEBDA \\
    \bottomrule
    \end{tabular}%
  \label{tab:tab1}%
    \caption{Summary of numerical statistics for the four pairs of networks ($m=1,2,3,4$). { Details of computation and comparison of statistical distributions are provided in the main text and relevant figures.}}
\end{sidewaystable}%

We select cases with $m=2$ and $m=3$ to visualize the curve of these statistics for the pair of networks in Fig.\ref{stat2}.
Although not shown, we note here that the case with $m=1$ is similar to the case $m=2$, while $m=4$ is  similar to $m=3$.

\begin{figure}
\begin{center}
\includegraphics[width=0.65\textwidth]{figs/fig4a.png}
\includegraphics[width=0.65\textwidth]{figs/fig4b.png}
\end{center}
\caption{Distribution of the statistics of networks with $m=2$ (upper four panels) and $m=3$ (lower four panels). Diameter is a straight histogram count, and the probability densities in the other panels are computed with a Gaussian kernel estimator. We see systematic bias and decreased variance for preferential attachment across all indicators --- exactly as one would expect. Most importantly, PA networks are smaller (diameter) typically more strongly dissassortative (see \cite{sJ10}). The ${\rm SF}(\gamma,m)${} networks are labelled in the figure captions as UniS --- for uniform random sampling.}
\label{stat2}
\end{figure}

From the curves of  various network statistics depicted above, we find that PA networks are atypical. In particular, PA networks have: (i)  more negative assortativity compared with corresponding ${\rm SF}(\gamma,m)${} networks; (ii) increasing clustering coefficient as $m$ increases; and, (iii)  in the case $m=3$ and $m=4$ the two clustering coefficient curves separate from each other. All these  discrepancies point to the fact that as $m$ increases, the PA network becomes more and more hub-centric, while the ${\rm SF}(\gamma,m)${} network remain highly uniform -- and the high degree nodes  remain evenly distributed throughout the ${\rm SF}(\gamma,m)${} networks.

We now introduce algorithms to modify the specific aspects of the network
which contribute to this hub-centric property. Note that the hub-centric
structure is a global property of the network, since only modification on
a small portion of PA networks (such as the case with the so-called
``rich-club'' phenomena \cite{xu4}) is not sufficient: when we only
make the modification described in \cite{xu4} to manipulate rich-club
connections the PA and ${\rm SF}(\gamma,m)${} network statistical distributions remain
disparate. Therefore, the modification scheme we propose in the following
applies across the entire network structure --- from super-rich nodes and
hubs, to the poorest nodes. To maintain conciseness and focus, we present
the brief idea of our modifications scheme here; for details see
Algorithm~\ref{alg:spread} of  \ref{modification}. In the
following paragraphs we outline these two modification schemes --- one to
remove what we call the hub-centric features of PA networks and a second
scheme to add these features to ${\rm SF}(\gamma,m)${} realisations. The aim of these
computations is to show that these modifications alone are sufficient to
align the corresponding distributions of network structural properties
(assortativity and so on).

The modification scheme for an initial PA network is the following. First, we cut the links within the group of rich nodes and between rich nodes and the group of poor nodes as far as possible while preserving the connectivity of the network. Then, we use the algorithm described in \cite{xu4} among the giant nodes (which we define as the nodes whose degree is even larger than the typical rich nodes), to reconstruct the structure of the rich nodes. That is, we obtain networks with minimal interconnection between hubs, minimal connection from hubs to low degree nodes, but a similar rich-club structure (the interconnection among the super hubs).  After setting the thresholds $\alpha_1=60\%$, $\alpha_2=5\%$ and $\alpha_3=0.5\%$ in Alg.~\ref{alg:spread} (\ref{modification}), we apply this modification to the case with $m=3$ and $m=4$. The result with $m=3 $ is shown in Fig.~\ref{mod1}, and for conciseness, the result with $m=4$, although similar, are omitted.

\begin{figure}
\begin{center}
\includegraphics[width=0.95\textwidth]{figs/fig5
\end{center}
\caption{(\textbf{Modification on PA}) The probability distribution function curve of the statistics for ${\rm SF}(\gamma,m)${} networks (labelled UniS, solid lines), PA networks (bold dotted gray lines), modified PA networks (termed PA\'\, with bold dotted black lines) and modified fitted PA networks (termed PA\'\,\'\,with dotted black lines), where $m=3$.}
\label{mod1}
\end{figure}

From Fig.~\ref{mod1}, the range of the {{probability distribution functions}} for the modified graphs overlap with the unmodified target distribution ${\rm SF}(\gamma,m)${} --- indicating that our modification provides a good fit to the distribution of PA networks, with the {{addition}} of hub-centric links. Moreover, we test the modified curves of the PA networks and the curves of ${\rm SF}(\gamma,m)${} networks. For this purpose, we add the following random factors. First, we choose $\alpha_1$ randomly from the set $[50\%,80\%]$, and then the result for $m=3 $ is also depicted in Fig.~\ref{mod1}. The results for $m=4$ are similar, but omitted for conciseness.

From the calculations of Fig.~\ref{mod1} we see that the modification of the PA scheme (removing the hub-centric properties of such networks) allows us to produce network statistics which have a distribution similar to the unmodified uniform sampling scheme (for the local clustering coefficient and assortativity \cite{sJ10}, at least) or bracketing the expected distribution for ${\rm SF}(\gamma,m)${} (for the maximum diameter). Although this bracketing --- by modifying the hub-centric nature of the PA network we go from smaller than ${\rm SF}(\gamma,m)${} to larger that ${\rm SF}(\gamma,m)${} --- does not immediately provide indistinguishable statistical distributions, it is clear that this is due only to the unselective manner in which we choose the threshold parameters for this algorithm. Changing these parameters effects a continuous change in these statistical distributions. This is sufficient to make our case that these hub-centric groups of nodes are what causes the PA networks to be atypical.  We do note, however, that the parametric changes we have explored are insufficient to modify the PA network  to reproduce the full distribution of observed global clustering coefficients for the ${\rm SF}(\gamma,m)${} network: there appears to be still further unexplored richness in the variety of these networks --- scale-free networks with close to zero clustering (very tree-like networks \cite{kJ13}).

\begin{figure}
\begin{center}
\includegraphics[width=0.95\textwidth]{figs/fig6}
\end{center}
\caption{(\textbf{Modification on ${\rm SF}(\gamma,m)${}}) The probability distribution function curve of the ${\rm SF}(\gamma,m)${} networks (labelled UniS with solid lines), PA networks (bold dotted gray lines) and modified ${\rm SF}(\gamma,m)${} networks (termed UniS' with dotted lines) in the case $m=3$}
\label{rev3}
\end{figure}

To further support our claim, we also present here a reverse modification
algorithm on ${\rm SF}(\gamma,m)${} networks to generate interconnection among hubs,
which would result in the adjusted curves of the statistics for ${\rm SF}(\gamma,m)${}
networks approaching that of the PA networks. Leaving the details to
\ref{modification} (Alg.~\ref{alg:concentrate}), we again
introduce the idea briefly here.  First, we delete the edges among the
poor nodes and relink these edges to rich nodes with probability
proportional to the degree of the rich nodes. Second, we use the
algorithm described in \cite{xu4} to make a club of super-rich
(``giant'') nodes (the definition is similar to the one in the first
modification) connected. In this way we are able generate the hub-centric
structure in networks. Setting $\alpha_1=80\%$, $\alpha_2=5\%$ and
$\alpha_3=0.5\%$ for Alg.~\ref{alg:spread}, the result for $m=3$ is shown
in Fig.~\ref{rev3}. The calculation for $m=4$ is similar and omitted for
conciseness.

From Fig.~\ref{rev3}, after generating hubs by artificially reconstructing the structure of networks, the numerical analysis of ${\rm SF}(\gamma,m)${} networks indicates statistical distributions of structural properties which approach that of PA networks --- supporting our claim that hub-centric structure in PA networks makes them atypical of random scale-free graphs.

By modifying the properties of a uniformly sampled network we are able to
produce distributions of maximum diameter and global clustering
coefficient very similar to that observed for the PA. That is, we add
hub-centric structure to  ${\rm SF}(\gamma,m)${} network and achieve statistical
features similar to that typical of PA networks. We do not achieve such
good agreement for assortativity and local clustering
coefficient. Nonetheless, taken together Fig.~\ref{mod1} and~\ref{rev3}
show that matching distributions can be achieved by {\em either} adding or
removing the hub-centric features in each of the four statistics that we
examine here. Hence, the features which we modify with the algorithms
described in \ref{modification} are exactly the properties of PA networks that
are atypical of the broader distribution of properties of scale free networks
when sampled by a maximum entropy process. The only caveat being that the
two algorithms are not exactly inverses of one another, and, in
particular, measures of global clustering hint at further unexplored
diversity in the global structure of typical scale free networks (as the
networks become progressively more tree like, even for $m>1$).

\subsection{Hub-centric structure in PA networks}

To conclude this section, we now collate the results of our analysis.  We
see that the hubs of a PA network ``hold'' that network together.  The
strength of interconnection among hubs (and connection with low degree
nodes) may be explained by preferential attachment itself. In that
growth model, once nodes and links are added, they will never be altered
again. This, we claim, causes inhomogeneity in PA networks. The largest
hubs will always (with probability approaching one) be the earliest nodes
in the growing network and these nodes are necessarily closely
interlinked. Conversely, the last nodes added to the network will have
minimal degree and yet these nodes will (with very high probability) be
directly connect to the hubs. In fact, it is clear that the last nodes
added (those with the lowest degree) are {\em most} likely to be connected
directly to the hubs  --- this is a consequence of the fact that the
  attachment bias favouring the original rich nodes increases with growth
  of the network --- loosely speaking, that the ``rich get richer''. These
are precisely the properties we alluded to earlier with our description of
``hub-centric'' networks, and, these are precisely the properties which
are adjusted by the modification algorithms~\ref{alg:spread} and~\ref{alg:concentrate}
of \ref{modification}.

The consequences of this hub-centric structure of PA networks is two-fold. {\em First}, since PA is an elegant and intuitive way to generate graphs, there may exist some generation procedures in the real world which are similar to preferential attachment, and thus we can use this claim to illustrate the hub-centric structure in such growth networks. This observation has potential for practical use --- for example, with hub detection in control of disease transmission, or,  to control of a network by manipulating the hubs.
{\em Second}, such a claim also indicates that there is systematic bias with PA networks. This bias will lead to difficulty when using PA networks to explain real world networks which do not result from such a constrained growth process --- even when the degree sequence of each network satisfies the power law distribution.

\section{Surrogate Networks: An application to particular putative scale-free systems}
\label{sec:applications}

In this section, we introduce a variant of the surrogate data test,
proposed for nonlinear time series analysis, to interpret real world
networks. We proceed by generating an ensemble of random networks, both
similar to the observed data-based network (in that they share the same
degree distribution) and yet, at the same time, random. By comparing the
properties of the real network with the corresponding distribution of
${\rm SF}(\gamma,m)${} networks we can determine whether the particular network is
typical. Our observations here are both a consequence of the previous
section and a motivation for the development of a network surrogate test.

Just like the surrogate data methods for time series, results will depend on the choice of test statistic. We will see that certain test statistics indicate clear and systematic deviation between the maximum entropy realisations and particular experimental systems. This can be interpreted in two ways. {\em First}, as a straight-forward hypothesis test. Rejection on the basis of any particular statistic indicates that the observed network is not typical of the maximum entropy class. However, we are also able to be more particular. Hence, {\em second}, we can apply this method to determine which particular features of the observed network are atypical. Hence we provide a more sensitive classification scheme dependent on individual features of the network. In what follows we will see that this is very often the case --- it is particular features of the networks that differ, rather than {\em all} aspects {\em simultaneously.}

We will start with a brief discussion of robustness. As we saw above, PA
networks are vulnerable to targeted attack, while ${\rm SF}(\gamma,m)${} networks don't
have such an evident ``Achilles' heel''. Recent research into the
structure of several important complex networks shows that, even if their
degree distribution could have an approximately power-law tail, the
networks in question are robust to targeted attack to some degree: the most
highly connected nodes do not necessarily represent an ``Achilles' heel''. In particular, recent results of modeling the router-level
Internet has shown that the core of that network is constructed from a
mesh of high bandwidth, low-connectivity routers, and \cite{jD05b}
concludes that although power-law degree distributions imply the presence
of high-degree vertices, they do not imply that such nodes form a
necessarily ``crucial hub''. A related line of research into the structure
of biological metabolic networks have shown that claims of SF structure
fail to capture the most essential biochemical as well as ``robust yet
fragile'' features of cellular metabolism and in many cases completely
misinterpret the relevant biology --- for example, \cite{rT05}
indicates the removal of high degree nodes leaves the biosynthetic
pathways fully intact. Hence, real-world scale-free networks do, indeed,
exhibit absence of the much-touted ``robustness'' properties. Our model
provides an explanation for that absence.

In the following, we will probe the limitation of the explanative power of PA networks via numerical statistics. We find that  many real complex networks appear more ``uniform'' under our surrogate test. Also, based on this result, we propose a simplistic classification of real world networks. Such a classification could provide us with a new insight into the interior structure of real networks, to explore which property do make the networks interesting and special.

Our surrogate tests are developed in this way: for a real world network
which satisfies the power-law distribution, we estimate its minimum vertex
degree and power law exponent by using the algorithm described in
\cite{aC09}. We then use the relation $\gamma=2+A/m$ in {Section
  $1$} to estimate its initial attractiveness, and then generate PA
networks with controlled $\gamma$ and $m$. We then substitute $\gamma$ and
$m$ into Algorithm~\ref{alg:MEp} to generate the ${\rm SF}(\gamma,m)${} networks. For
each set of parameter values, we generate $40$ networks for ${\rm SF}(\gamma,m)${} and PA
networks respectively.  In Fig.~\ref{box1}, we draw the boxplot of
numerical statistics of the real networks, as well as its corresponding PA
and ${\rm SF}(\gamma,m)${} networks.

\begin{figure}[!htp]
\begin{center}

\includegraphics[width=0.95\textwidth]{figs/fig7.pdf}
\end{center}
\caption{Network characteristics: (a) diameter; (b) global clustering coefficient; (c) local clustering coefficient; and, (d) assortativity. Boxplot analysis for collaboration and information networks, depicting  maximum, minimum, upper quantile ($75\%$), median and lower quantile ($25\%$) of the data.  $1.$ CS PhD collaboration\, $2.$ Erd\"os collaboration \, $3.$ a symmetrized snapshot of the structure of the Internet at the level of autonomous systems\, $4.$ a metabolic pathway network,  $5.$ the {\em S.cerevisiae} protein-protein interaction network\, $6.$ US Airport connection.}
\label{box1}
\end{figure}

\begin{figure}[!hp]
\begin{center}
\includegraphics[width=0.65\textwidth]{figs/fig8a.png}
\includegraphics[width=0.90\textwidth]{figs/fig8b.png}
\end{center}
\caption{The first panel shows the motif distribution and the lower panels the curves of statistics of the ${\rm SF}(\gamma,m)${} networks (solid lines, labelled UniS), PA networks (dotted lines) and  Erd\"os collaboration network (bold mixed lines). The motif rank in ascending order is shown below the graph.}
\label{erdosmotif}
\end{figure}

From Fig. \ref{box1} we can conclude that, in terms of interpretation of real networks, ${\rm SF}(\gamma,m)${} networks provides a significantly better model for  data from the Internet topology, and metabolic and protein interaction networks. For CS PhD collaboration and US airport connection networks, ${\rm SF}(\gamma,m)${} is slightly better than PA. However, for Erd\"os collaboration network, neither ${\rm SF}(\gamma,m)${} nor PA networks is good, although the statistics of PA networks are closer to that of the Erd\"os network than ${\rm SF}(\gamma,m)${} networks.

To examine the case of Erd\"os collaboration network more closely, we plot the motif rank and various other statistics of the Erd\"os collaboration networks in Fig. \ref{erdosmotif}.
This figure suggests that the hub-centric structure in the Erd\"os collaboration network is much more significant than PA allows: the Erd\"os collaboration network is a densely connected core along with loosely coupled radial branches reaching out from that core.  Erd\"os practiced what he preached -- he was a weaver of social networks and thus a builder of social capital. Moreover, this collaboration network is specifically Erd\"os-centric --- it is specifically focussed on that one unique hub and its connections. Hence, the nodes closer to Erd\"os benefit most and become the strongest hubs (after Erd\"os himself) in the resultant network.

By comparison, it would be helpful for us to identify which properties differ most markedly between the real network and a random surrogate sample.. The cases where ${\rm SF}(\gamma,m)${} provides better representations would suggest, for
instance, that the structure of protein interaction networks is more
uniform, because the cellular ecosystem necessitates such stableness; and,
in the case of the the Internet on the level of automatic system, the
system is balanced and distributed in a rather deliberate way. Both these
systems have been ``engineered'' for robustness.  Finally, we note that
neither the ${\rm SF}(\gamma,m)${} nor PA networks performs well in reproducing the
diameter of the biological networks. It is curious that biological
networks present abnormally large diameters --- we speculate that this
more spread-out tree-like structure is a consequence of the components of
biological systems having specific functional roles.

Based on our results, we also note (via analysis of robustness, motif and numerical statistics of a network, and comparison with its corresponding PA and ${\rm SF}(\gamma,m)${} networks) that one can classify the real network in the sense of hub-structure: if its property is closer to the PA networks, we can conclude it is hub-centric and typically the result of a growth process. Conversely, networks that are not PA are more likely to be purposefully widely distributed.

\section{Conclusion and discussion}
\label{conclude}

We have proposed a new algorithm which allows us to make maximum-entropy
random-samples of finite size graphs with the degree histogram being a
probabilistic realization of a specified degree distribution.  We focus
on  putative scale free networks as an illustrative example and use this method to generate random networks
with power-law tail degree distributions of arbitrary $\gamma$ (exponent)
and $m$ (minimum degree). This provides a simple alternative to the
various generative processes in the literature. However, our approach has
the benefit that it makes no particular biasing assumptions (such a
preferential attachment). To emphasize the need for this algorithm we
compare our results to distributions of networks obtained from the widely
applied PA scheme of \cite{aB99}.  While there are many generative algorithms for scale-free networks, PA is still the most widely used. Hence, PA provides an excellent model of
growth of a static network, however, we find that many real world networks
do not conform to this model.  This is not a new observation, however, it is still an important point to emphasise --- growth (and preferential attachment in particular) are not adequate models to explain all observed scale-free networks. We  propose particular structural modification algorithms to uncover exactly how PA is atypical of what one would expect of maximum entropy.

We found the high-degree nodes in PA networks are hub-centric, and that this hub-centric-ness has a greater influence on the overall structure than attributable only to the high-degree nodes themselves. In particular we attribute the co-called ``robust yet fragile'' feature found in many scale free networks to the distribution of hubs achieved via preferential attachment. Our observation sod the hub-centric-ness of preferential attachment scale free networks helps us assess the contribution of hubs in real world networks to the overall network structure. The surrogate network generation method we describe allows us to provide a robust numerical measure of this rich variation among scale-free complex networks.