\section{Introduction}

Random walks often provide the underlying
mesoscopic mechanism for transport phenomena
in physics, chemistry and biology \cite{MoSh84,MeKl00,KlSo11}.
A wide class of random walks give rise to normal diffusion,
where the mean-square displacement (MSD),
$\langle (\Delta r)^{2}(t)\rangle$,
grows linearly with time $t$ for long times.
In many important applications, however, the MSD behaves
like $\langle (\Delta r)^{2}(t)\rangle\propto t^{\gamma}$,
with $\gamma\neq 1$,
and the diffusion is anomalous \cite{MoSh84,MeKl00}.
Anomalous diffusion can be modelled by various classes of random
walks \cite{MeKl04}. We focus on the important class of
continuous time random walks (CTRWs) \cite{MoSh84,MeKl00}.
A specific feature of a CTRW is that
a walker waits for a random time
$\tau$ between any two successive jumps.
These waiting times are random independent
variables with a probability distribution
function (PDF) $\phi(\tau)$, and the tail of the PDF
determines if the transport is
diffusive ($\gamma=1$) or subdiffusive ($\gamma < 1$).
Heavy-tailed waiting time PDFs give rise to
subdiffusion.
Realistic models of the
waiting time PDF have been formulated for transport in
disordered materials with fractal and ramified
architecture, such as porous discrete media \cite{MaBaShBr03}
and
comb and dendritic polymers \cite{CaBe66,DoRoFr90,Fr05},
and for transport in crowded environments \cite{So12}.

A simple caricature of various types
of natural branched structures that
belong to the category of loopless graphs
is a comb model (see Fig.~\ref{fig:f1}).
The comb model was introduced to understand anomalous
transport in percolation clusters \cite{WhBa84,WeHa86,ArBa91}.
Now, comb-like models are widely employed
to describe various experimental
applications. These models have proven useful
to describe the transport
along spiny dendrites \cite{MeIo13,IoMe13},
percolation clusters with dangling bonds \cite{WeHa86},
diffusion of
drugs in the circulatory system \cite{MaRiMc08},
energy transfer in comb
polymers \cite{CaBe66,DoRoFr90} and
dendritic polymers \cite{Fr05}, diffusion
in porous materials \cite{ArKuBa11,StCo84,TaJa12},
the  influence of
vegetation architecture on the diffusion of
insects on plant surfaces
\cite{Ha02}, and many other interdisciplinary applications.
Random walks on comb structures provide
a geometrical explanation
of anomalous diffusion.

More general combs have been studied recently.
For example,
a numerical study of the encounter problem of
two walkers
in branched structures shows that the topological
heterogeneity of
the structure can play an important role \cite{AgBlCa14}.
Another example
is the occupation time statistics for random walkers
on combs
where the branches
can be regarded as independent complex structures,
namely fractal or other
ramified branches \cite{ReBa13}. Finally, we want to
mention studies
to understand the diffusion mechanism along
a variety of branched
systems, where scaling arguments, verified by
numerical simulations, have
been able to predict how the MSD
grows with time \cite{FoBuCeVu13}.

Diffusion on comb structures has also been
studied by
macroscopic approaches,
based on
Fokker-Planck equations \cite{ArBa91},
which have been applied
to describe diffusive properties
in discrete systems, such as
porous discrete media \cite{MaBaShBr03},
infiltration of diffusing particles
from one material into another \cite{KoBa10}, and
superdiffusion due
to the presence of inhomogeneous
convection flow \cite{BaIo04,IoBa05}.
Other macroscopic descriptions, based
on renormalizing the waiting
time PDF for jumps along
the backbone
to take into account the
transport along the branches
\cite{vdb89}, have been found
useful
to model
continuous-time-reaction-transport
processes \cite{CaMe05} and
human migrations along river networks \cite{CaFoMe06}.

Kahng and Redner provided a mesoscopic,
probabilistic description of random walks
on combs, by using the
successive decimation of the discrete-time Master
equation to obtain a
mesoscopic balance equations for the probability of
the walker to be at
a given node at a given time \cite{KaRe89}.
A mesoscopic description
is necessary for
 an accurate description of the transport
properties, such as
the diffusion coefficient or the mean visiting
time in a branch, in terms
of the microscopic parameters that
characterize the random walk
process.

Here we obtain transport
quantities within the framework of the CTRW
formalism. We assume that the
motion along the backbone
and the branches is non-Markovian and
that the motion
along the branches can be non-isotropic.
We reduce the dynamic effect of the branches
to a
waiting time PDF
for the motion along the backbone by using the
decimation method of Kahng and Redner.
The time spent by the walker between its entry
into the branches and its return to
the backbone for the first time is
treated as a contribution to the
effective waiting time at the node
where the branch crosses the backbone.

The paper is organized as follows.
In Sec.~\ref{sec:mdesc} we formulate
the mesoscopic description of the random walk
on the comb and reduce walker's motion
to an effective motion
along the backbone only with
a renormalized waiting time PDF for
the backbone nodes. Sec.~\ref{sec:stat} deals
with the MSD of the effective backbone motion,
derives the effective diffusion coefficient,
and establishes the conditions
for normal diffusion, anomalous
diffusion, and stochastic localization
(diffusion failure) \cite{DeHo00}
in terms of
the number of branch nodes and the degree
of bias of the motion along the branches.
We provide details of the numerical calculations
in Sec.~\ref{sec:numer} and summarize our
results and discuss their
implications in Sec.~\ref{sec:concl}.

\section{Mesoscopic description}\label{sec:mdesc}

The simplest comb model, shown in Fig. \ref{fig:f1},
is formed by a principal axis, called the
backbone, which is a one-dimensional
lattice with spacing $a$,  and identical branches
that
cross the backbone perpendicularly at each node.
\begin{figure}[htbp]
\includegraphics[width=0.7\hsize]{./figure1.eps}
\caption{Comb structure consisting of a backbone
and branches. Each point represents a
node where the walker may jump or wait a random time}
\label{fig:f1}
\end{figure}
The walker moves through the comb by
performing jumps between
nearest-neighbor nodes along
the backbone or along the branches.
We assume that
the walker performs isotropic jumps along the backbone,
but
the
jumps along the branches may be biased,
for example by an external field \cite{WhBa84}.

We derive the balance equation for the PDF $P(x,t)$
of finding
the walker at node $x$ on the backbone at time $t$.
When the walker arrives at a node,
it waits a random time $\tau$
before performing a new jump to the nearest node.
We assume that the comb is homogeneous,
and the waiting time PDF at any given node
is given by $\phi_{0}(\tau)$.
When the walker enters a branch, it spends some time
moving inside the branch
before returning to the backbone.
This sojourn time can be used to
determine an effective waiting time PDF
$\phi(\tau)$ for the walker's motion along
the backbone.
In other words,
the motion of the walker on the comb
can be reduced
to the effective
motion along a one-dimensional lattice,
corresponding to the backbone only.
This motion is non-Markovian and can be
described mesoscopically by
the Generalized Master Equation
(GME)
\begin{multline}
\frac{\partial P}{\partial t}=\int_{0}^{t}K(t-t')dt'\times\\
\left[\int_{-\infty}^{\infty} P(x-x',t')\Phi(x')dx'-P(x,t')\right],
\label{eq:me}
\end{multline}
where $K(t)$ is the memory kernel related to the waiting time PDF
via its Laplace transform,
$K(s)=s{\phi}(s)/[1-{\phi}(s)]$,
where $s$ is the Laplace variable.
The dispersal kernel $\Phi(x)$
represents the probability for the walker
of performing a jump of
length $x$. If the walker moves
isotropically between
nearest neighbors in a one-dimensional
lattice of spacing $a$,
the dispersal kernel reads
$\Phi(x)=\delta(x-a)/2+\delta(x+a)/2$.
We assume that the walker is initially
located at $x=0$, i.e., $P(x,0)=\delta_{x,0}$
with $x=ia$ and $i=0,\pm 1,\pm 2, \dotsc$, where
$\delta_{x,0}$ is the Kronecker delta.
Then the Laplace transform of
the GME for $x\neq 0$
reads
\begin{equation}
P(x,s)=\frac{\phi(s)}{2}\left[P(x-a,s)+P(x+a,s)\right].
\label{eq:me2}
\end{equation}

The mesoscopic balance equation for the walker
on the comb
being at node $x=ia$
of the backbone is
\begin{multline}\label{eq:me3}
P(x,s)=\frac{\phi_{0}(s)}{4}\left[P(x-a,s)+P(x+a,s)\right] \\
+(1-q)\phi_{0}(s)\left[P(y=a,s)+P(y=-a,s)\right].
\end{multline}
Here $P(x,s)$, $P(x-a,s)$, and $P(x+a,s)$
is short-hand for
$P(x,y=0,s)$, $P(x-a,y=0,s)$, and $P(x+a,y=0,s)$,
and $P(y=a,s)$ and $P(y=-a,s)$ stands for
$P(x,y=a,s)$ and $P(x,y=-a,s)$.
The term $\phi_{0}(s)\left[P(x-a,s)+P(x+a,s)\right]/4$
corresponds to
the contribution of
the walker arriving at node $x=ia$ from the left
or from the right with probability 1/4
after waiting a random time $\tau$ with PDF
$\phi_{0}(\tau)$ at nodes $x+a$ or $x-a$.
As shown in Fig.~\ref{fig:f2},
the walker located at
the $i$th node of
the backbone may jump to the right, left, up
or down with probability
1/4. We assume that the walker moves forward
(away from the backbone)
along the
branches with probability $q$ and back to the
backbone with probability $1-q$.
The term
\begin{equation}
(1-q)\phi_{0}(s)\left[P(y=a,s)+P(y=-a,s)\right]
\end{equation}
in (\ref{eq:me3}) corresponds the contribution of
the walker arriving
at the backbone node $x$
from the first node of the upper or
lower branch after waiting
there a random time
$\tau$ with PDF
$\phi_{0}(\tau)$.

Consider the motion along the upper branches.
The lower branch
dynamics is the same due to the symmetry of the comb.
The mesoscopic balance
equation for the first node of the upper branches reads
\begin{equation}
\label{eq:me4}
P(y=a,s)=\frac{\phi_{0}(s)}{4}P(x,s)+\phi_{0}(s)(1-q)P(y=2a,s).
\end{equation}
The first term $\phi_{0}(s)P(x,s)/4$ corresponds to
the
contribution of the walker
arriving from the backbone,
while $\phi_{0}(s)(1-q)P(y=2a,s)$ is
the contribution of the walker jumping
from the upper node $y=2a$
to $y=a$ with probability $1-q$
after waiting  a random time
$\tau$ with PDF
$\phi_{0}(\tau)$.
Analogously, we have for the lower branches
\begin{equation}
\label{eq:me4b}
P(y=-a,s)=\frac{\phi_{0}(s)}{4}P(x,s)+\phi_{0}(s)(1-q)P(y=-2a,s).
\end{equation}
Generalizing (\ref{eq:me4}) to any node of the branches
located between  $2a\leq y\leq(N-2)a$,
we obtain the balance equation
for the upper branches
\begin{equation}
\label{eq:me5}
P(y,s)=\phi_{0}(s)\left[qP(y-a,s)+(1-q)P(y+a,s)\right].
\end{equation}

To determine the Laplace transform $\phi(s)$
of the effective backbone node waiting
time PDF, we need to determine
$P(y=a,s)$ and $P(y=-a,s)$
in \eqref{eq:me3}
in terms of $P(x,t)$, so that
\eqref{eq:me3} can be cast in the form
of \eqref{eq:me2}. Given \eqref{eq:me4}
and \eqref{eq:me4b}, this goal can be achieved
if $P(y=2a,s)$ and $P(y=-2a,s)$
can be related
to $P(y=a,s)$ and $P(y=-a,s)$.
We proceed as follows.
The solution of (\ref{eq:me5}) reads
\begin{equation}
\label{eq:solme5}
P(y,s)=A_{1}\lambda_{+}^{y/a}+A_{2}\lambda_{-}^{y/a},
\end{equation}
where
\begin{equation}
\lambda_{\pm}=\frac{1\pm\sqrt{1-4q(1-q)\phi_{0}^{2}(s)}}
{2(1-q)\phi_{0}(s)}.
\end{equation}
To find expressions for
the quantities $A_{1}$ and $A_{2}$,
whose dependence on $x$ and $s$ is not displayed,
we apply
(\ref{eq:solme5}) to the node $y=2a$:
\begin{equation}\label{eq:sol2me5}
P(y=2a,s)=A_{1}\lambda_{+}^{2}+A_{2}\lambda_{-}^{2}.
\end{equation}
On the other hand, setting $y=2a$ in (\ref{eq:me5}),
we find
\begin{multline}
P(y=2a,s)=\phi_{0}(s)\left[qP(y=a,s)\right.\\
+\left.(1-q)\phi_{0}(s)P(y=3a,s)\right],
\end{multline}
or
\begin{multline}
P(y=2a,s)-\phi_{0}(s)qP(y=a,s)=\\
\phi_{0}(s)(1-q)\phi_{0}(s)P(y=3a,s).
\end{multline}
Setting $y=3a$
in (\ref{eq:solme5}) we obtain
\begin{multline}
\label{eq:sol3me5}
P(y=2a,s)-q\phi_{0}(s)P(y=a,s)\\
=\phi_{0}(s)(1-q)\left[A_{1}\lambda_{+}^{3}+A_{2}\lambda_{-}^{3}\right].
\end{multline}
Solving the system of equations (\ref{eq:sol2me5}) and (\ref{eq:sol3me5})
for the quantities $A_{1}$ and $A_{2}$, we obtain
\begin{multline}
\label{eq:a1}
A_{1}=\frac{P(y=2a,s)-q\phi_{0}(s)P(y=a,s)}
{\lambda_{+}^{2}\left(\lambda_{+}-\lambda_{-}\right)\phi_{0}(s)(1-q)}\\
-\frac{\lambda_{-}P(y=2a,s)}{\lambda_{+}^{2}\left(\lambda_{+}-\lambda_{-}\right)},
\end{multline}
\begin{multline}
\label{eq:a2}
A_{2}=\frac{-P(y=2a,s)+q\phi_{0}(s)P(y=a,s)}
{\lambda_{-}^{2}\left(\lambda_{+}-\lambda_{-}\right)\phi_{0}(s)(1-q)}\\
+\frac{\lambda_{+}P(y=2a,s)}{\lambda_{-}^{2}\left(\lambda_{+}-\lambda_{-}\right)}.
\end{multline}

A special situation occurs at the end of the branches,
where we have
to impose reflecting boundary conditions, i.e.,
\begin{equation}
\label{eq:me6}
P(y=Na,s)=q\phi_{0}(s)P(y=(N-1)a,s).
\end{equation}
The node at $y=(N-1)a$ also needs a special balance equation
(see
Fig.~\ref{fig:f2}),
\begin{figure}[htbp]
\includegraphics[width=0.5\hsize]{./figure2.eps}
\caption{Schematic representation of the possible jumps
of a walker with the corresponding probabilities.}
\label{fig:f2}
\end{figure}
\begin{multline}
\label{{eq:me7}}
P(y=(N-1)a,s)=q\phi_{0}(s)P(y=(N-2)a,s)\\
 +\phi_{0}(s)P(y=Na,s).
\end{multline}

Substituting $y=(N-2)a$ into (\ref{eq:me5})
and considering (\ref{eq:me6}),
we can write
\begin{equation}
\label{eq:p2p3}
P(y=(N-2)a,s)=h(\phi_{0}(s))P(y=(N-3)a,s),
\end{equation}
where
\begin{equation}
\label{eq:h}
h(\phi_{0}(s))=\frac{q\phi_{0}(s)\left[1-q\phi_{0}^{2}(s)\right]}
{1+q(q-2)\phi_{0}^{2}(s)}.
\end{equation}
Substituting
the solutions from
(\ref{eq:solme5}), (\ref{eq:a1}), and (\ref{eq:a2})
into (\ref{eq:p2p3}),
we find
\begin{equation}
\label{eq:p2p1}
P(y=2a,s)=G(q,\phi_{0}(s))P(y=a,s),
\end{equation}
where
\begin{align}
\label{eq:gh}
G(q,\phi_{0}(s))&=\frac{2q\phi_{0}(s)}{1+
\dfrac{1+H(q,\phi_{0}(s))}{1-H(q,\phi_{0}(s))}
\sqrt{1-4q(1-q)\phi_{0}^{2}(s)}},\\
H(q,\phi_{0}(s))&=\left(\frac{\lambda_{-}}{\lambda_{+}}\right)^{N-5}
\frac{\lambda_{-}-h(\phi_{0}(s))}{\lambda_{+}-h(\phi_{0}(s))}.
\end{align}
For the lower branch we obtain in
a similar manner,
\begin{equation}
\label{eq:pm2m1}
P(y=-2a,s)=G(q,\phi_{0}(s))P(y=-a,s).
\end{equation}
We have achieved our goal of
expressing
$P(y=2a,s)$ and $P(y=-2a,s)$
in terms of $P(y=a,s)$ and $P(y=-a,s)$.
Substituting  (\ref{eq:p2p1}) and (\ref{eq:pm2m1})
into (\ref{eq:me4})
and (\ref{eq:me4b}) and
using the resulting expressions in (\ref{eq:me3}),
we obtain an equation of the
form (\ref{eq:me2}) with
\begin{equation}
\label{eq:wt}
\phi(s)=\frac{\phi_{0}(s)}{2-
\dfrac{(1-q)\phi_{0}^{2}(s)}{1-(1-q)\phi_{0}(s)G(q,\phi_{0}(s))}}.
\end{equation}
The Laplace inversion of (\ref{eq:wt}) yields
$\phi (\tau)$, which incorporates the dynamics
along the branches and can be understood as the
effective waiting time PDF for a walker
moving along the backbone only.

\section{Statistical properties}\label{sec:stat}
\subsection{$N$ finite}

If the local waiting time PDF $\phi_{0}(\tau)$ has
finite moments, its Laplace transform
reads \cite{MeKl00}, $\phi_{0}(s)\simeq1-s\bar{t}$,
in the large time limit $s\rightarrow0$,
where $\bar{t}$ is the local mean waiting
time at each node. Taking the limit
$s\rightarrow0$ in (\ref{eq:wt}),
we obtain the waiting time PDF
for the effective backbone
dynamics,
\begin{equation}\label{eq:phi}
\phi(s)\simeq(1+s\left\langle t\right\rangle )^{-1}.
\end{equation}
The mean waiting time $\left\langle t\right\rangle$
is given by
\begin{equation}
\label{eq:mwt}
\left\langle t\right\rangle =
\frac{\bar{t}}{2q-1}\left[2(1-q)^{1-N}q^{N}+4q-3\right].
\end{equation}
\begin{figure}[htbp]
\includegraphics[width=0.9\hsize]{./figure3.eps}
\caption{Dimensionless
mean waiting time of the
effective backbone dynamics.}
\label{fig:f3}
\end{figure}
In Fig.~\ref{fig:f3},
we plot the effective mean waiting at a node
of the backbone dynamics
versus $N$ and $q$.
It shows that the
mean waiting time $\left\langle t\right\rangle$
is a monotonically increasing function of both $q$ and $N$.
If the random walk inside the branches is isotropic,
$q=1/2$,
one obtains by the L'Hopital's rule from (\ref{eq:mwt})
\begin{equation}
\lim_{q\rightarrow1/2}\left\langle t\right\rangle
=\left(1+2N\right)\bar{t}.
\end{equation}
To determine the diffusion coefficient $D$
for diffusion through the comb,
we first calculate the MSD.
Performing the Fourier-Laplace transform
on (\ref{eq:me}), we obtain
\begin{equation}\label{eq:me-FL}
P(k,s)=\frac{1-\phi(s)}{s[1-\Phi(k)\phi(s)]}.
\end{equation}
The MSD in Laplace
space reads (see, e.g., \cite{MeKl00})
\begin{equation}
\label{eq:msd}
\left\langle x^{2}(s)\right\rangle =
-\lim_{k\rightarrow0}\frac{d^{2}P(k,s)}{dk^{2}}.
\end{equation}
As mentioned in Sec.~\ref{sec:mdesc},
we assume that the motion on the
backbone is unbiased and that the walker
only jumps to nearest neighbors.
This implies that the kernel $\Phi(x)$
is given $\Phi(x)=\delta(x-a)/2+\delta(x+a)/2$,
and we obtain from
(\ref{eq:msd}),
\begin{equation}
\label{eq:msd2}
\left\langle x^{2}(s)\right\rangle
=\frac{a^{2}}{s\left[\phi(s)^{-1}-1\right]}\,
\end{equation}
in the large time limit.
If
the waiting time PDF
$\phi(t)$ possesses a finite first moment,
(\ref{eq:phi}) implies that the
MSD along the backbone corresponds
to normal diffusion
$\left\langle x^{2}(t)\right\rangle =2Dt$.
The diffusion coefficient is given by
\begin{equation}\label{eq:d}
 D=\frac{a^{2}}{2\left\langle t\right\rangle}
 =\frac{a^{2}}{2\bar{t}}\frac{2q-1}{2(1-q)^{1-N}q^{N}+4q-3}.
\end{equation}
In Fig.~\ref{fig:f4}, we compare the
results provided by (\ref{eq:d})
with numerical simulations.
\begin{figure}[htbp]
\includegraphics[width=1.\hsize]{./figure4.eps}
\caption{Plot of the diffusion coefficient through the comb
for $N=2$, $N=3$, and $N=6$ versus $q$.
Solid curves correspond to exact analytical
results given by (\ref{eq:d}).
Numerical simulations results are depicted with symbols.}
\label{fig:f4}
\end{figure}

As Fig.~\ref{fig:f3} demonstrates,
$\left\langle t\right\rangle$
increases monotonically with $N$
for $q<1/2$ and saturates at
$(4q-3)/(2q-1)$
for $N\to\infty$. Consequently, the
mean waiting time $\left\langle t\right\rangle$
is finite for $N\to\infty$;
the overall
diffusion along the backbone
is normal.
However, for $q\geq1/2$, the mean waiting time
$\left\langle t\right\rangle$
increases without bound as $N$ increases,
and anomalous transport
is expected for $N\to\infty$.
\begin{figure}[htbp]
\includegraphics[width=1.\hsize]{./figure5.eps}
\caption{MSD over $2D$ for $N$ fixed and three different
values of $q$: 0.1, 0.25, 05. 0.75}
\label{fig:f5}
\end{figure}
In Fig.~\ref{fig:f5}, we plot the MSD scaled
by the diffusion coefficient.
It illustrates the result given by (\ref{eq:d}) for the MSD.
The transport is diffusive for finite $N$, regardless $q$ and
the specific form of $\phi_{0}(\tau)$, as
long as it has finite moments.

We consider now a an effective
waiting time PDF with
the large-time limit
$\phi_{0}(\tau)\sim \tau^{-1-\gamma}$,
with Laplace transform
$\phi_{0}(s)\simeq 1-(s\tau_{0})^{\gamma}$
and $0<\gamma<1$, which
does not possess finite moments.
Here $\tau_{0}$ is a parameter with units of time.
In this case, the waiting time PDF
for the backbone dynamics is obtained by
simply replacing $s\bar{t}$ with $(s\tau_{0})^{\gamma}$,
i.e, $\phi(s)\simeq[1+\left(s\tau_{0}\right)^{\gamma}
\left\langle t\right\rangle /\tau_{0}]^{-1}$.
Substituting this result into (\ref{eq:msd2}),
we find
\begin{equation}
\left\langle x^{2}(t)\right\rangle
=\frac{a^{2}\tau_{0}}{\left\langle t\right\rangle }
\frac{(t/\tau_{0})^{\gamma}}{\Gamma(1+\gamma)},
\end{equation}
for large $t$, where $\left\langle t\right\rangle $
is given by (\ref{eq:mwt}), with $\tau_{0}$ instead
of
$\bar{t}$.
If the waiting time PDF $\phi_{0}(\tau)$
at each node of the comb has a
power-law tail,
then the overall transport along
the backbone is anomalous.

\subsection{$N\to\infty$}

If the number of nodes of the branches
goes to infinity,
the mean time spent by the walker visiting
a branch increases monotonically,
see \eqref{eq:mwt}.
However, this does not always
results in anomalous transport
along the overall structure as
we show below.

For $N\to\infty$,
the quotient $(\lambda_{-}/\lambda_{+})^{N}\to 0$
and also $H\to 0$. We obtain from
(\ref{eq:gh}),
\begin{equation}
G(q,\phi_{0}(s))=\frac{2q\phi_{0}(s)}
{1+\sqrt{1-4q(1-q)\phi_{0}^{2}(s)}}\equiv
\frac{2q\phi_{0}(s)}{1+g(q)},
\end{equation}
where we define
$g(q)\equiv\sqrt{1-4q(1-q)\phi_{0}^{2}(s)}$
for  convenience.
Equation (\ref{eq:wt}) reduces to
\begin{equation}
\label{eq:fis}
\phi(s)=\frac{\phi_{0}(s)
\left[1+g(q)-2q(1-q)\phi_{0}^{2}(s)\right]}
{2-(1+3q-4q^2)\phi_{0}^{2}(s)+
[2-(1-q)\phi_{0}^{2}(s)]g(q)}.
\end{equation}
We take the limit $s \to 0$ and
consider first the case where
$\phi_{0}(\tau)$ has finite moments.
Then $\phi_{0}(s)\simeq 1-s\bar{t}$,
as $s\to 0$. The square root $g(q)$ in
(\ref{eq:fis}) reads
\begin{equation}
g(q)\simeq
\begin{cases}
1-2q-\dfrac{4q(1-q)}{2q-1}s\bar{t}, & q<1/2,\\[2ex]
\sqrt{2s\bar{t}}-\dfrac{\sqrt{2}}{4}(s\bar{t})^{3/2}, & q=1/2,\\[2ex]
-1+2q+\dfrac{4q(1-q)}{2q-1}s\bar{t}, & q>1/2,
\end{cases}
\end{equation}
and (\ref{eq:fis}) implies that the waiting time PDF
is given by
\begin{equation}
\label{eq:fis2}
\phi(s)\simeq
\begin{cases}
\left(1+\dfrac{4q-3}{2q-1}s\bar{t}\right)^{-1}, & q<1/2,\\[2ex]
\left(1+\sqrt{2s\bar{t}}\right)^{-1}, & q=1/2,\\[2ex]
\left(\dfrac{3q-1}{q}+\dfrac{4q^{2}-3q+1}{(2q-1)q}s\bar{t}\right)^{-1}, & q>1/2.
\end{cases}
\end{equation}
Substituting (\ref{eq:fis2}) into (\ref{eq:msd2}),
we find for large $t$,
\begin{equation}
\label{eq:msdi}
\left\langle x^{2}(t)\right\rangle =
\begin{cases}
a^{2}\dfrac{2q-1}{4q-3}\dfrac{t}{\bar{t}}, & q<1/2,\\[2ex]
a^{2}\sqrt{\dfrac{2t}{\pi\bar{t}}}, & q=1/2,\\[2ex]
a^{2}\dfrac{q}{2q-1}\left(1-e^{-\alpha t}\right),& q>1/2,
\end{cases}
\end{equation}
where the rate of saturation is
\begin{equation}
\alpha = \frac{(2q-1)^{2}}{(4q^{2}-3q+1)\bar{t}}.
\end{equation}
In Fig~\ref{fig:f6} we compare these
results with numerical simulations
for $N=10^{3}$.
For $q=1/2$, we obtain the
well known result of subdiffusive transport
with the MSD $\sim \sqrt{t}$.
However, for $q\neq 1/2$, the side branches experience
advection, and the transport is remarkably different.
Namely, for $q>1/2$ the advection is
away from the backbone
along
the branches, $y\rightarrow\pm\infty$.
The walker is effectively trapped
inside the branches,
and
stochastic localization (diffusion failure)
occurs,
$\left\langle x^{2}(\infty)\right\rangle
< \infty$,
\cite{DeHo00}.
For $q<1/2$, the advection
is towards the backbone.
It enhances the backbone dynamics
and normal diffusion takes place.
\begin{figure}[htbp]
\includegraphics[width=1.\hsize]{./figure6.eps}
\caption{MSD for three values of $q$,
displaying three different behaviors.
Solid curves correspond to the results
given by (\ref{eq:msdi}).
Symbols are the results of numerical
simulations with $N=10^3$.}
\label{fig:f6}
\end{figure}

Consider now the case where the
local waiting time PDF is $\phi_{0}(\tau)
\sim \tau^{-1-\gamma}$, i.e.,
$\phi_{0}(s)\simeq1-(s\tau_{0})^{\gamma}$
with $0<\gamma<1$, as $s\to 0$.
The MSD in this case
can be obtained straightforwardly
by replacing $s\bar{t}$ with $(s\tau_{0})^{\gamma}$
in (\ref{eq:fis2}). For large times it reads
\begin{equation}
\label{eq:msd-gamma}
\left\langle x^{2}(t)\right\rangle =
\begin{cases}
\dfrac{a^{2}}{\Gamma(1+\gamma)}\dfrac{2q-1}{4q-3}\left(\dfrac{t}{\tau_{0}}\right)^{\gamma}, & q<1/2,\\[2ex]
\dfrac{a^{2}}{\sqrt{2}\Gamma(1+\gamma/2)}\left(\dfrac{t}{\tau_{0}}\right)^{\gamma/2}, & q=1/2,\\[2ex]
\dfrac{a^{2}q(2q-1)}{\tau_{0}^{\gamma}(4q^{2}-3q+1)}\mu(t/\tau_{0}), & q>1/2,
\end{cases}
\end{equation}
where
\begin{equation}
\mu(t/\tau_{0})=(t/\tau_{0})^{\gamma}
E_{\gamma,\gamma+1}\left[-\left(\frac{t}{\tau_{0}}\right)^{\gamma}
\frac{(2q-1)^{2}}{4q^{2}-3q+1}\right]
\end{equation}
is expressed in terms
of the the generalized
Mittag-Leffler function $E_{\alpha,\beta}(z)$.
We use the following property of integration
of the Mittag-Leffler function \cite{Po99},
\begin{equation}
\int_0^tE_{\alpha,\beta}\left(bz^{\alpha}\right)z^{\beta-1}dz=t^{\beta}
E_{\alpha,\beta+1}\left(bt^{\alpha}\right).
\end{equation}
Subdiffusion in the branches results in
backbone subdiffusion for $q\leq 1/2$.
For advection away from the backbone,
$q>1/2$, we again
find stochastic localization.
For $t/\tau_{0}\gg 1$,
$E_{\alpha,\beta}(-at^{\alpha})
\sim t^{-\alpha}/\Gamma(\beta-\alpha)$
\cite{Ba53}, and
consequently
$\mu(t/\tau_{0})$ approaches
a finite value as $t\to\infty$.

\section{Conclusion}\label{sec:concl}

We have developed a mesoscopic equation
for a random walk on a regular comb structure
given by (\ref{eq:me2}) and (\ref{eq:wt}).
The random walk along the branches consists
of, possibly biased, jumps to the
nearest node, while waiting
at each node for a random time $\tau$
distributed according to the PDF $\phi_0(\tau)$
before proceeding with the next jump.
The overall dynamics along the branches
has been reduced to an
effective waiting time PDF,
given by (\ref{eq:wt}),
for motion solely
along the backbone. We have obtained
statistical properties, such as the
effective mean
waiting time, $\langle t\rangle$ for
the backbone nodes,
and the diffusion coefficient, $D$, of
the overall structure for the case where
the number of nodes $N$ of the branches
is finite or infinite.
If $N$ is finite and $\phi_0(\tau)$
has finite moments, both $\langle t\rangle$ and
$D$ are derived exactly in terms of the bias
probability $q$, the number
of nodes $N$ on the branch,
and the mean waiting time
probability at each node.
In this case the
transport is always
normal diffusion.
If $\phi_0(\tau)\sim \tau^{-1-\gamma}$
for large time, it does not
posses finite moments and the
MSD of the random walker
behaves like $t^{\gamma}$.
If $N$ is infinite,
the value of $q$ is decisive.
If $\phi_0(\tau)$ has finite moments,
the diffusion
regime is normal if $q<1/2$, while the
MSD behaves like $t^{1/2}$ for $q=1/2$.
If $q>1/2$,
the MSD approaches a constant finite
value
for large time, corresponding
to stochastic localization (diffusion
failure).
If $\phi_0(\tau)\sim \tau^{-1-\gamma}$
for large time, the MSD behaves
like $t^{\gamma}$
for $q<1/2$ and like $t^{\gamma/2}$ for $q=1/2$
Again, stochastic localization
occurs for $q>1/2$.
In summary, if the bias probability of moving
away from the backbone is $q>1/2$, then
stochastic localization occurs,
regardless of the other
characteristic parameters related to the
random walk on the branches.

\section*{Acknowledgments}

A.I. would like to thank the Universitat Aut\`onoma
de Barcelona for hospitality and financial support,
as well as the support by the
Israel Science Foundation (ISF-1028). VM and DC have
been supported by the Ministerio de Ciencia e
Innovaci\'on under Grant No. FIS2012-32334.
VM also thanks the Isaac Newton Institute for
Mathematical Sciences, Cambridge, for support
and hospitality during the CGP programme where
part of this work was undertaken.


\section{Introduction and preliminaries}
 In this paper we give new criteria for the verification of higher order convex orders. These criteria can be used to prove the Hermite-Hadamard type inequalities for higher order convex functions.

Let $f \colon [a,b]\to\mbox{$\mathbb{R}$}$ be a convex function. The following double inequality
\begin{equation}
f\left(\frac{a+b}2\right)\leqslant \frac1{b-a}\int_a^b f(x)\,dx \leqslant \frac{f(a)+f(b)}2 \label{eq:oli0}
\end{equation}
is known as the Hermite-Hadamard inequality for convex functions (see  \cite{DraPe2000} for many generalizations and applications of \eqref{eq:oli0}).

In many papers (see, for example, \cite {Florea, Rajba12f, Rajba14, Rajba12g, Szostok2014a, TSzostok2014, OlbrysSzostok2014}) are studied the Hermite-Hadamard type inequalities based on the convex stochastic ordering properties. In the paper \cite{Rajba12f}, to get a simple proof of some known Hermite-Hadamard type inequalities as well as to obtaining new Hermite-Hadamard type inequalities, is used the Ohlin lemma  on sufficient conditions for convex stochastic ordering. Recently, the Ohlin lemma is also used  to study the inequalities of the Hermite-Hadamard type in \cite {Rajba14, Rajba12g, Szostok2014a, TSzostok2014, OlbrysSzostok2014}. In the papers \cite {Szostok2014a, TSzostok2014, OlbrysSzostok2014}, furthermore, to examine the Hermite-Hadamard type inequalities is used the Levin-Ste\v{c}kin theorem \cite {LevinSteckin1960} (see also \cite {NicPer06}), which gives necessary and sufficient conditions for the stochastic convex ordering.

Let us recall some basic notions and results on the stochastic convex order (see, for example, \cite {DenLefSha98}). As usual, $F_X$ denotes the distribution function of a random variable $X$ and $\mu_X$ is the distribution corresponding to $X$. For real valued random variables $X,Y$ with a finite expectation, we say that $X$ is dominated by $Y$ in \textit{convex ordering} sense if
 $$\mbox{$\mathbb{E}$} f(X) \leqslant \mbox{$\mathbb{E}$} f(Y)$$
for all convex functions $f \colon \mbox{$\mathbb{R}$} \to \mbox{$\mathbb{R}$}$ (for which the expectations exist). In that case we write $X \leqslant_{\textit{cx}} Y$, or $\mu_X \leqslant_{\textit{cx}} \mu_Y$.

In the following Ohlin lemma \cite{Ohlin69} are given sufficient conditions for convex stochastic ordering.
\begin{lemma}[\cite{Ohlin69}]\label{lemma:3}
Let $X,Y$ be two random variables such that $\mbox{$\mathbb{E}$} X=\mbox{$\mathbb{E}$} Y$. If the distributions functions $F_X, F_Y$ cross exactly one time, i.e., for some $x_0$ holds
\begin{equation}\label{eq:ls1.2}
F_X(x) \leqslant F_Y(x) \textit{ if } x < x_0 \textit{ and } F_X(x) \geqslant F_Y(x) \textit{ if } x > x_0,
\end{equation}
then
\begin{equation}\label{eq:levin1}
\mbox{$\mathbb{E}$} f(X) \leqslant \mbox{$\mathbb{E}$} f(Y)
\end{equation}
for all convex functions $f \colon \mbox{$\mathbb{R}$} \to \mbox{$\mathbb{R}$}$.
\end{lemma}
\begin{remark}
As noticed Szostok in \cite{Szostok2014a}, if the measures $\mu_X $, $\mu_Y$ corresponding to $X$ and $Y$, respectively, are concentrated on the interval $[a,b]$ then, in fact, inequality \eqref{eq:levin1} is satisfied for all continuous convex functions $f \colon [a,b] \to \mbox{$\mathbb{R}$}$ .
\end{remark}
\begin{remark}
The inequality \eqref{eq:oli0} may be easily proved with the use of the Ohlin lemma (see\cite{Rajba12f}). Indeed, let $X$, $Y$, $Z$ be three random variables with the distributions $\mu_X = \delta_{(a+b)/2}$, $\mu_Y$ which is equally distributed in $[a,b]$ and $\mu_Z = \frac{1}{2}(\delta_a+\delta_b)$, respectively. Then it is easy to see that the pairs $(X,Y)$ and $(Y,Z)$ satisfy the assumptions of the Ohlin lemma, and using \eqref{eq:levin1}, we obtain \eqref{eq:oli0}.
\end{remark}
As we can see, the Ohlin lemma is a strong tool, however, it is worth noticing that in the case of some inequalities, the distribution functions cross more than once. Therefore a simple application of the Ohlin lemma is impossible and some additional idea is needed.

In the papers \cite{TSzostok2014, OlbrysSzostok2014}, the authors used the Levin-Ste\v{c}kin theorem  \cite{LevinSteckin1960} (see also \cite{NicPer06}, Theorem 4.2.7).
\begin{theorem}[\cite{LevinSteckin1960}]\label{th:1.3}
Let $a,b\in\mbox{$\mathbb{R}$}$, $a<b$ and let $F_1,F_2 \colon [a,b]\to\mbox{$\mathbb{R}$}$ be functions with bounded variation such that $F_1(a)=F_2(a)$. Then,
in order that
\begin{equation*}
\int_a^bf(x)dF_1(x)\leqslant\int_a^bf(x)dF_2(x),
\end{equation*}
for all continuous convex functions $f \colon [a,b]\to\mbox{$\mathbb{R}$},$
it is necessary and sufficient that $F_1$ and $F_2$ verify the following three conditions:
\begin{eqnarray}
F_1(b)& = & F_2(b)\label{eq:ls1.4},\\
\int_a^b F_1(x)dx & = & \int_a^b F_2(x)dx, \label{eq:ls1.6}\\
  \int_a^x F_1(t)dt & \leqslant & \int_a^x F_2(t)dt\quad for\: all\:x\in(a,b).\label{eq:ls1.5}
\end{eqnarray}
\end{theorem}
\begin{remark}
Observe, that if the measures $\mu_X$, $\mu_Y$, corresponding to the random variables $X$, $Y$, respectively, occurring in the Ohlin lemma, are concentrated on the interval $[a,b]$, the Ohlin lemma is an easy consequence of Theorem \ref{th:1.3}. Indeed, $\mu_X$ and $\mu_Y$ are probabilistic measures, thus we have $F_X(a)=F_Y(a)=0$ and $F_X(b)=F_Y(b)=1$. Moreover, $\mbox{$\mathbb{E}$} X=\mbox{$\mathbb{E}$} Y$ yields \eqref{eq:ls1.6} and from the inequalities \eqref{eq:ls1.2} we obtain \eqref{eq:ls1.5}.

\end{remark}

Szostok \cite{TSzostok2014} used Theorem \ref{th:1.3} to make an observation, which is more general than Ohlin's lemma and concerns the situation when the functions $F_1$ and $F_2$ have more crossing points than one. First we need the following definitions.

Define the number of sign changes of a function $\varphi \colon \mbox{$\mathbb{R}$} \to \mbox{$\mathbb{R}$}$ by
$$
S^-(\varphi) = \sup \{S^-[\varphi(x_1),\varphi(x_2),\ldots,\varphi(x_k)]\colon x_1<x_2<\ldots x_k\in\mbox{$\mathbb{R}$},\: k\in\mbox{$\mathbb{N}$}\},
$$
where $S^-[y_1,y_2,\ldots,y_k]$ denotes the number of sign changes in the sequence $y_1,y_2,\ldots,y_k$ (zero terms are being discarded). Two real functions $\varphi_1,\varphi_2$ are said to have $n$ crossing points (or cross each other $n$-times) if $S^-(\varphi_1-\varphi_2)=n$. Let $a=x_0<x_1< \ldots<x_n<x_{n+1}=b$. We say that the functions  $\varphi_1,\varphi_2$ crosses $n$-times at the points $x_1,x_2,\ldots,, x_n$ (or that $x_1,x_2,\ldots,, x_n$ are the points of sign changes of $\varphi_1-\varphi_2$) if $S^-(\varphi_1-\varphi_2)=n$ and there exist $a<\xi_1<x_1< \ldots<\xi_n<x_n<\xi_{n+1}<b$ such that $S^-[\xi_1,\xi_2,\ldots,\xi_{n+1}]=n$.

Then, Lemma 2 given in \cite{TSzostok2014} can be rewritten in the following form.
\begin{lemma}[\cite{TSzostok2014}]\label{lemma:1.5}
Let $a,b\in\mbox{$\mathbb{R}$}$, $a<b$ and let $F_1,F_2 \colon (a,b)\to\mbox{$\mathbb{R}$}$ be functions with bounded variation such that $F_1(a)=F_2(a)$, $F_1(b)=F_2(b)$, $F=F_2-F_1$, $\int_a^b F(x)dx=0$. Let $a<x_1< \ldots<x_m<b$ be the points of sign changes of $F$ and $F(t)\geq 0$   for $t\in(a,x_1)$.
\begin{itemize}
\item If $m$ is even then
the inequality
\begin{equation}\label{eq:ord1}
\int_a^bf(x)dF_1(x)\leqslant\int_a^bf(x)dF_2(x)
\end{equation}
is not satisfied by all continuous convex functions $f \colon [a,b]\to\mbox{$\mathbb{R}$}.$

\item If $m$ is odd, define $A_i$ ($i=0,1,\ldots,m$, $x_0=a$, $x_{m+1}=b$)
$$
A_i = \int_{x_i}^{x_{}i+1} |F(x)|dx.
$$
Then the inequality \eqref{eq:ord1} is satisfied for all continuous convex functions $f \colon [a,b]\to\mbox{$\mathbb{R}$},$  if and only if the following inequalities hold true:
\begin{equation}\label{eq:a1}
\begin{split}
A_0 & \geqslant  A_1, \\
A_0+A_2 & \geqslant  A_1+A_3 ,\\
& \vdots \\
A_0+A_2+ \ldots+A_{m-3} & \geqslant  A_1+A_3+ \ldots+A_{m-2}.
\end{split}
\end{equation}
\end{itemize}
\end{lemma}
\begin{remark}\label{rmk:a2}
Let
\begin{equation*}
H(x)=\int_a^x F(t)dt.
\end{equation*}
Then the inequalities  \eqref{eq:a1}
are equivalent to the following inequalities
\begin{equation*}
H(x_2)\geqslant  0,\: H(x_4)\geqslant  0,\:H(x_6)\geqslant  0, \:\ldots, \:H(x_{m-1})\geqslant 0.
\end{equation*}
\end{remark}

Now we are going to study Hermite-Hadamard type inequalities for higher-order convex functions. Many results on higher order generalizations of the Hermite-Hadamard type inequality one can found, among others, in \cite{Bes08, BesPal02, BesPal03, BesPal04, BesPal10, DraPe2000, Rajba12f, Rajba12g}. In recent papers \cite{Rajba12f, Rajba12g} the theorem of M. Denuit, C.Lefevre and  M. Shaked \cite{DenLefSha98} on sufficient conditions for $s$-convex ordering was used to prove Hermite-Hadamard type inequalities for higher-order convex functions.

Let us review some notations. The convexity of $n$-th order (or $n$-convexity) was defined in terms of divided differences by Popoviciu \cite{Popoviciu1934}, however, we will not state it have. Instead we list some properties of $n$-th order convexity which are equivalent to Popoviciu's definition (see \cite{Kuczma1985}).

\begin{proposition}\label{prop:16}
A function $f \colon (a,b)\to\mbox{$\mathbb{R}$}$ is $n$-convex on $(a,b)$ $(n \geqslant 1)$ if and only if its derivative $f^{(n-1)}$ exists and is convex on $(a,b)$ (with the convention $f^{(0)}(x) = f(x)$).
\end{proposition}

\begin{proposition}\label{prop:17}
Assume that $f \colon [a,b] \to\mbox{$\mathbb{R}$}$ is $(n+1)$-times differentiable on $(a,b)$ and continuous on $[a,b]$ ($n \geqslant 1$). Then $f$ is $n$-convex if and only if $f^{(n+1)}(x)\geqslant 0$, $x \in (a,b)$.
\end{proposition}

For real valued random variables $X,Y$ and any integer $s \geqslant 2$ we say that $X$ is dominated by $Y$ in $s$-\textit{convex ordering} sense if $\mbox{$\mathbb{E}$} f(X) \leqslant \mbox{$\mathbb{E}$} f(Y)$ for all $(s-1)$-convex functions $f \colon \mbox{$\mathbb{R}$} \to \mbox{$\mathbb{R}$}$, for which the expectations exist (\cite{DenLefSha98}). In that case we write $X \leqslant_{s-cx} Y$, or $\mu_X \leqslant_{s-cx} \mu_Y$, or $F_X \leqslant_{s-cx} F_Y$. Then the order $\leqslant_{2-cx}$ is just the usual convex order $\leqslant_{cx}$.

A very useful criterion for the verification of the $s$-convex order is given by Denuit, Lef\`{e}vre and Shaked in \cite{DenLefSha98}.
\begin{proposition}[\cite{DenLefSha98}]\label{prop:18}
Let $X$ and $Y$ be two random variables such that $\mbox{$\mathbb{E}$} (X^j-Y^j)=0$, $j=1,2,\ldots,s-1$ ($s \geqslant 2$). If $S^-(F_X-F_Y)=s-1$ and the last sign of $F_X-F_Y$ is positive, then $X \leqslant_{s-cx} Y$.
\end{proposition}
Proposition \ref{prop:18} can be rewritten in the following form.
\begin{proposition}[\cite{DenLefSha98}]\label{prop:19}
Let $X$ and $Y$ be two random variables such that
\begin{equation*}
\mbox{$\mathbb{E}$} (X^j-Y^j)=0,\quad j=1,2,\ldots,s \:(s \geqslant 1).
\end{equation*}
If the distribution functions $F_X$ and $F_Y$ cross exactly $s$-times at points $x_1<x_2< \ldots <x_s$ and
\begin{equation*}
(-1)^{s+1}\left(F_Y(x)-F_X(x)\right)\geqslant 0\quad for\: all\: x\leqslant x_1
\end{equation*}
then
\begin{equation}\label{eq:ord5}
\mbox{$\mathbb{E}$} f(X) \leqslant \mbox{$\mathbb{E}$} f(Y)
\end{equation}
for all $s$-convex functions $f \colon \mbox{$\mathbb{R}$} \to \mbox{$\mathbb{R}$}$.
\end{proposition}
\begin{remark}
Observe, that if the measures $\mu_X$, $\mu_Y$, corresponding to the random variables $X$, $Y$, respectively, occurring in Proposition \ref{prop:19}, are concentrated on the interval $[a,b]$, then in fact inequality \eqref{eq:ord5} is satisfied for all continuous $s$-convex functions $f \colon [a,b] \to \mbox{$\mathbb{R}$}$.
\end{remark}

Proposition \ref{prop:19} is a counterpart of the Ohlin lemma concerning convex ordering.  This proposition  gives sufficient conditions for $s$-convex ordering, and is very useful for the verification of higher order convex orders,  however, it is worth noticing that in the case of some inequalities, the distribution functions cross more than $s$-times. Therefore a simple application of this proposition is impossible and some additional idea is needed.

In this paper we give a theorem on necessary and sufficient conditions for higher order convex stochastic ordering, which can be useful in the study of Hermite-Hadamard type inequalities for higher order convex functions, and in particular inequalities between the quadrature operators. This theorem is a counterpart of the Levin-Ste\v{c}kin theorem \cite{LevinSteckin1960} concerning convex stochastic ordering. The necessary and sufficient conditions, which we give in this paper, concern higher order stochastic convex ordering of signed measures, and are generalization of results given by Denuit, Lef\`{e}vre and Shaked in  \cite{DenLefSha98}. Moreover, our criteria can be easier to checking of higher order convex orders, than those given in \cite{DenLefSha98}.
\section{Main results}
Let  $F_1,F_2 \colon \mbox{$\mathbb{R}$}\to\mbox{$\mathbb{R}$}$ be two functions with bounded variation and $\mu_1$, $\mu_2$ be the signed measures corresponding to $F_1$, $F_2 $, respectively. We say that $F_1$ is dominated by $F_2 $ in $(n+1)$-\textit{convex ordering sense} $(n\geqslant 1)$ if
\begin{equation*}
\int_{-\infty}^{\infty}f(x)dF_1(x)\leqslant\int_{-\infty}^{\infty}f(x)dF_2(x),
\end{equation*}
for all $n$-convex functions $f \colon \mbox{$\mathbb{R}$}\to\mbox{$\mathbb{R}$}$.
In that case we write $F_1 \leqslant_{(n+1)-cx} F_2$, or $\mu_1 \leqslant_{(n+1)-cx} \mu_2$.

In the following theorem we give necessary and sufficient conditions for $(n+1)$-convex ordering of two functions with bounded variation.
\begin{theorem}\label{th:2.1}
Let $a,b\in\mbox{$\mathbb{R}$}$, $a<b$, $n\in \mbox{$\mathbb{N}$}$ and let $F_1,F_2 \colon [a,b]\to\mbox{$\mathbb{R}$}$ be two functions with bounded variation such that $F_1(a)=F_2(a)$. Then,
in order that
\begin{equation*}
\int_a^bf(x)dF_1(x)\leqslant\int_a^bf(x)dF_2(x),
\end{equation*}
for all continuous $n$-convex functions $f \colon [a,b]\to\mbox{$\mathbb{R}$},$
it is necessary and sufficient that $F_1$ and $F_2$ verify the following  conditions:
\begin{equation}\label{eq:ls2.2}
F_1(b)= F_2(b),
\end{equation}
\begin{equation}\label{eq:ls2.3}
\int_a^b F_1(x)dx = \int_a^b F_2(x)dx,
\end{equation}
\begin{multline}\label{eq:ls2.4}
\int_a^b \int_a^{x_{k-1}} \ldots \int_a^{x_1} F_1(t)dtdx_1\ldots dx_{k-1}=\\
 \int_a^b \int_a^{x_{k-1}} \ldots \int_a^{x_1} F_2(t)dtdx_1\ldots dx_{k-1},\quad
 for\ \  k=2,\ldots ,n,
 \end{multline}
  \begin{multline}\label{eq:ls2.5}
(-1)^{n+1}\int_a^x \int_a^{x_{n-1}} \ldots \int_a^{x_1} F_1(t)dtdx_1\ldots dx_{n-1}\leqslant\\
 (-1)^{n+1}\int_a^x \int_a^{x_{n-1}} \ldots \int_a^{x_1} F_2(t)dtdx_1\ldots dx_{n-1}, \ \ for \ any \ x\in(a,b).
 \end{multline}
\end{theorem}

First we prove the following lemma.
\begin{lemma}\label{lemma:2.2}
Let $F \colon [a,b]\to\mbox{$\mathbb{R}$}$ be a function with bounded variation. Let $f \colon [a,b]\to\mbox{$\mathbb{R}$}$ be an $n$-convex function of the class $C^{n+1}$ on $(a,b)$. Then
\begin{equation}\label{eq:2.6}
\int_a^bf(x)dF(x)=\biggl[F(x)f(x)\biggl]_{x=a}^{x=b}-\int_a^bF(x)f'(x)dx,
\end{equation}
\begin{equation}\label{eq:2.7}
\int_a^bf(x)dF(x)=\biggl[F(x)f(x)\biggl]_{x=a}^{x=b}-\biggl[\int_a^xF(t)dtf'(x)\biggl]_{x=a}^{x=b}+\int_a^b\int_a^xF(t)dtf''(x)dx,
\end{equation}
\begin{multline}\label{eq:2.8}
\int_a^bf(x)dF(x)=\biggl[F(x)f(x)\biggl]_{x=a}^{x=b}-\biggl[\int_a^xF(t)dtf'(x)\biggl]_{x=a}^{x=b}+\biggl[\int_a^x\int_a^{x_1}F(t)dtdx_1f''(x)\biggl]_{x=a}^{x=b}+\ldots \\
\ldots +\biggl[(-1)^k\int_a^x\int_a^{x_{k-1}}\ldots  \int_a^{x_1}F(t)dtdx_1\ldots dx_{k-1}f^{(k)}(x)\biggl]_{x=a}^{x=b}+\\
+(-1)^{k+1}\int_a^b\int_a^x\int_a^{x_{k-1}}\ldots  \int_a^{x_1}F(t)dtdx_1\ldots dx_{k-1}f^{(k+1)}(x)dx,\quad for \:k=2,\ldots,n.
 \end{multline}
\end{lemma}
\begin{proof}
The proof is by induction. Integrating by parts and using the equalities $F(x)=\left( \int_a^xF(t)dt \right)'$ and $\int_a^xF(t)dt=\left(\int_a^x\int_a^{x_1}F(t)dtdx_1 \right)'$, we obtain immediately \eqref{eq:2.6}, \eqref{eq:2.7} and \eqref{eq:2.8} for $k=2$.

Put
$$
I_k(x)=\int_a^x\int_a^{x_{k}}\ldots  \int_a^{x_1}F(t)dtdx_1\ldots dx_{k},\quad for \:x\in(a,b), \:k=1,2,\ldots,n.
$$
Then we have
\begin{equation}\label{eq:2.9}
I_{k-1}(x)=\left(I_k(x)\right)',\quad for \:x\in(a,b), \:k=1,2,\ldots,n.
\end{equation}
Assume that \eqref{eq:2.8} holds for some $k=2,\ldots,n-1$. Integrating by parts and using \eqref{eq:2.9}, we obtain that the last summand in \eqref{eq:2.8} can be rewritten in the form
\begin{multline*}
(-1)^{k+1}\int_a^b\int_a^x\int_a^{x_{k-1}}\ldots  \int_a^{x_1}F(t)dtdx_1\ldots dx_{k-1}f^{(k+1)}(x)dx=\\=(-1)^{k+1}\int_a^bI_{k-1}(x)f^{(k+1)}(x)dx
=\biggl[I_k(x)f^{(k+1)}(x)\biggl]_{x=a}^{x=b}+(-1)^{k+2}\int_a^bI_{k}(x)f^{(k+2)}(x)dx,
 \end{multline*}
which implies that \eqref{eq:2.8} holds for $k+1$. Thus \eqref{eq:2.8} holds for all $k=2,\ldots,n$. The lemma is proved.
\end{proof}
The proof of Theorem \ref{th:2.1} is an immediate consequence of the following lemma.
\begin{lemma}\label{lemma:2.3}
Let $F \colon [a,b]\to\mbox{$\mathbb{R}$}$ be a function with bounded variation such that $F(a)=0$. Then in order that
\begin{equation}\label{eq:2.10}
\int_a^bf(x)dF(x)\geqslant 0,
\end{equation}
for all continuous $n$-convex functions $f \colon [a,b]\to\mbox{$\mathbb{R}$},$
it is necessary and sufficient that $F$ satisfies the following  conditions:
\begin{align}
F(b)&=0, \label{eq:2.11}\\
\int_a^b F(x)dx & =0,\label{eq:2.12}\\
\int_a^b \int_a^{x_{k-1}} \ldots \int_a^{x_1} F(t)dtdx_1\ldots dx_{k-1}&=0,\quad
 for\ \  k=2,\ldots ,n,\label{eq:2.13}\\
 (-1)^{n+1}\int_a^x \int_a^{x_{n-1}} \ldots \int_a^{x_1} F(t)dtdx_1\ldots dx_{n-1}&\geqslant 0,\quad for \ any \ x\in(a,b). \label{eq:2.14}
 \end{align}
\end{lemma}
\begin{proof}
Via an approximation argument we may restrict to the case when $f$ is of the class $C^{n+1}((a,b))$.

We now prove the sufficiency. By Lemma \ref{lemma:2.2}, using \eqref{eq:2.7} and \eqref{eq:2.8} with $k=n$, and taking into account \eqref{eq:2.11}-\eqref{eq:2.13} we get
\begin{equation}\label{eq:2.15}
\int_a^bf(x)dF(x)=(-1)^{n+1}\int_a^b\int_a^x \int_a^{x_{n-1}} \ldots \int_a^{x_1} F(t)dtdx_1\ldots dx_{n-1}f^{(n+1)}(x)dx.
\end{equation}
Then, by \eqref{eq:2.14} and Proposition \ref{prop:17}, we obtain \eqref{eq:2.10}.

We now prove the necessity. The necessity of \eqref{eq:2.11} follows by checking our statement for $f=1$ and $f=-1$.

 The necessity of \eqref{eq:2.12} follows by checking our statement for $f(x)=x$ and $f(x)=-x$ and by using \eqref{eq:2.11}, \eqref{eq:2.6}.

 The necessity of \eqref{eq:2.13} we prove by induction on $k$. The necessity of \eqref{eq:2.13} for $k=2$ follows by checking our statement for $f(x)=x^2$ and $f(x)=-x^2$, using \eqref{eq:2.7} and taking into account \eqref{eq:2.11}, \eqref{eq:2.12}. Assume, that the equality
 \begin{equation}\label{eq:2.16}
 \int_a^b \int_a^{x_{l-1}} \ldots \int_a^{x_1} F(t)dtdx_1\ldots dx_{l-1}=0
 \end{equation}
 holds , for some $k=2,\ldots ,n-1$ and all $l=2,\ldots ,k$. Then we check our statement for $f(x)=x^{k+1}$ and $f(x)=-x^{k+1}$. Using \eqref{eq:2.8} and taking into account \eqref{eq:2.11}, \eqref{eq:2.12} and \eqref{eq:2.16} for $l=2,\ldots ,k$, we obtain \eqref{eq:2.16} for $l=k+1$. Consequently, we obtain that \eqref{eq:2.13} is satisfied for all $k=2,\ldots ,n$.

 By \eqref{eq:2.8} with $k=n$ and taking into account \eqref{eq:2.11}-\eqref{eq:2.13}, we obtain that \eqref{eq:2.15} holds. Then, for the necessity of \eqref{eq:2.14}, notice that $(-1)^{n+1}\int_a^x \int_a^{x_{n-1}} \ldots \int_a^{x_1} F(t)dtdx_1\ldots dx_{n-1}<0$ for some $x\in(a,b)$, yields an interval $I$ around $x$  on which this expression is still negative. Choosing $f$ such that $f^{(n+1)}=0$ outside $I$, the equality \eqref{eq:2.15} leads to a contradiction. Thus \eqref{eq:2.14} is satisfied. The lemma is proved.
 \end{proof}
 \begin{corollary}\label{cor:2.4}
Let $\mu_1$, $\mu_2$ be two signed measures on $\mbox{$\mathcal{B}$}(\mbox{$\mathbb{R}$})$, which are concentrated on $(a,b)$, and such that $\int_a^b\lvert x \rvert ^n \mu_i(dx)<\infty$, $i=1,2$. Then in order that
\begin{equation*}
\int_a^bf(x)d\mu_1(x)\leqslant\int_a^bf(x)d\mu_2(x),
\end{equation*}
for continuous $n$-convex functions $f \colon [a,b]\to\mbox{$\mathbb{R}$}$, it is necessary and sufficient that  $\mu_1$, $\mu_2$ verify the following conditions:
\begin{align}
\mu_1\left((a,b)\right)&=\mu_1\left((a,b)\right),\label{eq:2.17}\\
\int_a^b x  ^k \mu_1(dx)&=\int_a^b x  ^k \mu_2(dx),\quad for \:k=1,\ldots ,n,\label{eq:2.18}\\
\int_a^b \bigl(t-x\bigr)  ^n _+\mu_1(dt)&=\int_a^b \bigl(t-x\bigr) ^n _+\mu_2(dt),\quad for \:x\in (a,b),\label{eq:2.19}
\end{align}
where $y^n_+=\Bigl\lbrace max \lbrace y,0\rbrace \Bigr\rbrace ^n$, $y\in \mbox{$\mathbb{R}$}$.
\end{corollary}
\begin{proof}
Let $F_1$, $F_2$ be the distribution functions corresponding to $\mu_1$, $\mu_2$, respectively. Then $\mu_i(dt)=dF_i(t)$, $i=1,2$. Since $\mu_1$ and $\mu_2$ are concentrated on $(a,b)$, we have $F_1(a)=F_2(a)$. That \eqref{eq:ls2.2} and \eqref{eq:2.17} are equivalent is obvious. Put $F=F_2-F_1$. By \eqref{eq:2.6} with $f(x)=x$, and taking into account \eqref{eq:ls2.2}, it follows that the conditions \eqref{eq:ls2.3} and \eqref{eq:2.18} for $k=1$ are equivalent. The equivalence of \eqref{eq:ls2.4} and \eqref{eq:2.18} for $k=2,\ldots ,n$, can be proved, using \eqref{eq:2.7} and \eqref{eq:2.18}, by induction on $k$. We omit the proof.

Next, by reversing the order of integration in \eqref{eq:ls2.5}, we obtain
\begin{align*}
&(-1)^{n+1}\int_a^x \int_a^{x_{n-1}} \ldots \int_a^{x_1} (F_2(t)-F_1(t))dtdx_1\ldots dx_{n-1}\\
&=(-1)^{n+1}\int_a^x \int_a^{x_{n-1}} \ldots \int_a^{x_1} F(t)dtdx_1\ldots dx_{n-1}\\
&=(-1)^{n+1}\int_a^x \dfrac{(x-t)^n}{n!}dF(t)=(-1)^{n+1}(-1)^{n}\int_a^x \dfrac{(t-x)^n}{n!}dF(t)=-\int_a^x \dfrac{(t-x)^n}{n!}dF(t)\\
&=\int_a^b \dfrac{(t-x)^n}{n!}dF(t)-\int_a^x \dfrac{(t-x)^n}{n!}dF(t)=\int_x^b \dfrac{(t-x)^n}{n!}dF(t)=\int_a^b \dfrac{(t-x)^n_+}{n!}dF(t),
\end{align*}
which implies the equivalence of \eqref{eq:ls2.5} and \eqref{eq:2.19}. The corollary is proved.
\end{proof}

Note that Theorem \ref{th:2.1} can be rewritten in the following form.
\begin{theorem}\label{th:2.7a}
Let $F_1,F_2 \colon [a,b]\to\mbox{$\mathbb{R}$}$ be two functions with bounded variation such that $F_1(a)=F_2(a)$. Let
\begin{align*}
H_0(t_0)&=F_2(t_0)-F_1(t_0), \quad for\: t_0\in [a,b],\\
H_{k}(t_k)&=\int_a^{t_{k-1}}H_{k-1}(t_{k-1})dt_{k-1},\quad for\: t_{k}\in [a,b], \  k=1,2,\ldots ,n.
\end{align*}
Then,
in order that
\begin{equation*}
\int_a^bf(x)dF_1(x)\leqslant\int_a^bf(x)dF_2(x),
\end{equation*}
for all continuous $n$-convex functions $f \colon [a,b]\to\mbox{$\mathbb{R}$},$
it is necessary and sufficient that the following  conditions are satisfied:
\begin{align*}
H_k(b)&=0, \quad for\: k=0,1,2,\ldots ,n,\\
(-1)^{n+1}H_{n}(x)&\geqslant 0 ,\quad for\:all\: x\in (a,b).
\end{align*}
\end{theorem}
\begin{remark}\label{rmk:2.8a}
The functions $H_1, \ldots,H_n$, that appear in Theorem \ref{th:2.7a} can be obtained from the following formulas
	\begin{equation}\label{eq:2.8c}
H_{n}(x)=(-1)^{n+1}\int_a^b \frac{(t-x)_+^n}{n!}d (F_2(t)-F_1(t)),
\end{equation}
\begin{equation}\label{eq:2.8d}
H_{k-1}(x)=H_{k}^{\,'}(x), \quad k=2,3,\ldots ,n.
\end{equation}
\end{remark}

Note that the function $(-1)^{n+1}H_{n-1}$ that appears in Theorem \ref{th:2.7a} play a role similar to the role of the function $F=F_2-F_1$ in Lemma \ref{lemma:1.5}. Consequently, from Theorem \ref{th:2.7a}, Lemma \ref{lemma:1.5} and Remarks \ref{rmk:a2}, \ref{rmk:2.8a} we obtain immediately the following useful criterion for the verification of higher order convex ordering.
\begin{corollary}\label{cor:2.7b}
Let $F_1,F_2 \colon [a,b]\to\mbox{$\mathbb{R}$}$ be functions with bounded variation such that $F_1(a)=F_2(a)$, $F_1(b)=F_2(b)$ and $H_{k}(b)=0$ $(k=1,2,\ldots ,n)$, where $H_k(x)$ $(k=1,2,\ldots ,n)$ are given by \eqref{eq:2.8c} and \eqref{eq:2.8d}. Let $a<x_1< \ldots<x_m<b$ be the points of sign changes of the function $H_{n-1}$ and let $(-1)^{n+1}H_{n-1}(x)\geqslant 0$   for $x\in(a,x_1)$.
\begin{itemize}
\item If $m$ is even then the inequality
 \begin{equation}\label{eq:2.7c}
\int_a^bf(x)dF_1(x)\leqslant\int_a^bf(x)dF_2(x),
\end{equation}
is not satisfied by all continuous $n$-convex functions $f \colon [a,b]\to\mbox{$\mathbb{R}$}$.

\item If $m$ is odd, then the inequality \eqref{eq:2.7c} is satisfied for all continuous $n$-convex functions $f \colon [a,b]\to\mbox{$\mathbb{R}$}$ if and only if
\begin{equation}\label{eq:2.7d}
(-1)^{n+1}H_n(x_2)\geq 0, (-1)^{n+1}H_n(x_4)\geq 0,\ldots,
 (-1)^{n+1}H_n(x_{m-1})\geq 0.
\end{equation}
\end{itemize}
\end{corollary}
In \cite{DenLefSha98} can be found the following necessary and sufficient conditions for the verification of the $(s+1)$-convex order.
\begin{proposition}[\cite{DenLefSha98}]\label{prop:2.5}
If $X$ and $Y$ are two real valued random variables such that $\mbox{$\mathbb{E}$} \lvert X \rvert^{s}<\infty$ and  $\mbox{$\mathbb{E}$} \lvert Y\rvert^{s}<\infty$ then
\begin{equation}\label{eq:2.5c}
\mbox{$\mathbb{E}$} f(X) \leqslant \mbox{$\mathbb{E}$} f(Y)
\end{equation}
for all continuous $s$-convex functions $f \colon \mbox{$\mathbb{R}$} \to \mbox{$\mathbb{R}$}$
if and only if
\begin{align}
\mbox{$\mathbb{E}$} X^k&=\mbox{$\mathbb{E}$} Y^k, \quad for \:k=1,2,\ldots,s,\label{eq:2.5a} \\
\mbox{$\mathbb{E}$} (X-t)^{s}_+&\leqslant \mbox{$\mathbb{E}$} (Y-t)^{s}_+, \quad for \:all\: t \in \mbox{$\mathbb{R}$}\label{eq:2.5b}.
\end{align}
\end{proposition}
\begin{remark}
The inequality \eqref{eq:2.5b} coincides with \eqref{eq:2.5c} for the spline function $f(x)=(x-t)^{s}_+$. Moreover, it is well known that $s$-convex function has the integral representation, such that the spline functions are the generating functions (see \cite{Rajba2011}).
\end{remark}
\begin{remark}
Note, that if the measures $\mu_X$, $\mu_Y$, corresponding to the random variables $X$, $Y$, respectively, occurring in Proposition \ref{prop:2.5}, are concentrated on some interval $[a,b]$, then this proposition is an easy consequence of Corollary \ref{cor:2.4}.
\end{remark}
\begin{remark}
Necessary and sufficient conditions for the verification of the $(s+1)$-convex order, which are given in Proposition \ref{prop:2.5}, can be difficult to checking. In this paper, we give conditions that can be easier for verification of higher order convex orders. Moreover, our conditions concern not only probabilistic measures but also signed measures.
\end{remark}

In the numerical analysis are studied some inequalities, which are connected with quadrature operators. These inequalities, called extremalities, are a particular case of the Hermite-Hadamard type inequalities. Many extremalities are known in the numerical analysis (cf. \cite{Bes08}, \cite{BraSch81}, \cite{BraPet} and the references therein).
The numerical analysts prove them using the suitable differentiability assumptions. As proved W\k{a}sowicz in the papers \cite{SzWas07b}, \cite{SzWas08}, \cite{SzWas10}, for convex functions of higher order some extremalities can be obtained without assumptions of this kind, using only the higher order convexity itself. The support-type properties play here the crucial role. As we show in \cite{Rajba12f, Rajba12g}, some extremalities can be proved using a probabilistic characterization.The extremalities, which we study are known, however our method using the Ohlin lemma \cite{Ohlin69} and the Denuit-Lef\`{e}vre-Shaked theorem \cite{DenLefSha98} on sufficient conditions for the convex stochastic ordering seems to be quite easy. It is worth noting that, these theorems do not apply to proving some extremalities (see \cite{Rajba12f, Rajba12g}). In these cases can be useful results given in this paper.

For a function $f:[-1,1]\to\mbox{$\mathbb{R}$}$ we consider two operators
\begin{eqnarray*}
 C(f)&:=&\tfrac{1}{3}\Bigl(
                      f\bigl(-\tfrac{\sqrt{2}}{2}\bigr)+f(0)+
					  f\bigl(\tfrac{\sqrt{2}}{2}\bigr)
                     \Bigr),\\\\
										\Lob{4}(f)&:=&\tfrac{1}{12}\Bigl(f(-1)+f(1)\Bigr)
              +\tfrac{5}{12}\Bigl(f\bigl(-\tfrac{\sqrt{5}}{5}\bigr)
              +f\bigl(\tfrac{\sqrt{5}}{5}\bigr)\Bigr),
\end{eqnarray*}
connected with Chebyshev and Lobatto quadratures, respectively.	W\k{a}sowicz \cite{SzWas07b}, \cite{SzWas08b} proved that
\begin{equation}\label{eq:w3}
 C(f)\leqslant \Lob{4}(f), \quad if \: f\:is \:3-convex.
\end{equation}
The proof given in \cite{SzWas07b} is rather complicated. This was done using computer software. In \cite{SzWas08b} can be found a new proof, based on the spline approximation of convex functions of higher order.
Using Corollary \ref{cor:2.7b} we give a new much simpler proof of \eqref{eq:w3}.

Since for the random variables $X$ and $Y$ with the distributions
\begin{eqnarray*}
\mu_X & = &\tfrac{1}{3}\Bigl( \delta_{-\frac{\sqrt{2}}{2}} + \delta_0+\delta_\frac{\sqrt{2}}{2}\Bigr), \\ \\
\mu_{Y} & = & \tfrac{1}{12}\Bigl(\delta_{-1}+\delta_1\Bigr)+\tfrac{5}{12}\Bigl(\delta_{-\frac{\sqrt{5}}{5}} + \delta_\frac{\sqrt{5}}{5}\Bigr),
\end{eqnarray*}
respectively, we have
$$
C(f) = E[f(X)],\quad \Lob{4}(f) = E[f(Y)],
$$
it follows that the inequality \eqref{eq:w3} can be rewritten in terms of higher order convex orderings
\begin{equation}\label{eq:w4}
 X \leqslant_{4-cx} Y.
\end{equation}
It is worth noting, that Proposition \ref{prop:2.5} of  Denuit, Lef\`{e}vre and Shaked does not apply to proving \eqref{eq:w4}, because of the distribution functions   $F_X$ and $F_Y$ cross exactly 5-times. We prove the inequality \eqref{eq:w4} by using Corollary \ref{cor:2.7b}.

We have $F_1=F_X$, $F_2=F_Y$, $H_0=F=F_Y-F_X$. By \eqref{eq:2.8c} and \eqref{eq:2.8d}, we obtain
\begin{eqnarray*}
H_3(x)&=&\tfrac{1}{72}\left\{ \left( -1-x \right)^3_+ +\left( 1-x \right)^3_+
+5\left[\left( -\tfrac{\sqrt{5}}{5}-x \right)^3_+ +\left( \tfrac{\sqrt{5}}{5}-x \right)^3_+\right] \right.\\ \\
\quad &\,& \left.-4\left[\left( -1-x \right)^3_+ +\left( -\tfrac{\sqrt{2}}{2}-x \right)^3_+ +\left( -x \right)^3_+ +\left( \tfrac{\sqrt{2}}{2}-x \right)^3_+ \right] \right\},
\end{eqnarray*}
\begin{eqnarray*}
H_2(x)&=&\tfrac{1}{24}\left\{ -\left( -1-x \right)^2_+ -\left( 1-x \right)^2_+
-5\left[\left( -\tfrac{\sqrt{5}}{5}-x \right)^2_+ +\left( \tfrac{\sqrt{5}}{5}-x \right)^2_+\right] \right.\\ \\
\quad &\,& \left.+4\left[\left( -1-x \right)^2_+ +\left( -\tfrac{\sqrt{2}}{2}-x \right)^2_+ +\left( -x \right)^2_+ +\left( \tfrac{\sqrt{2}}{2}-x \right)^2_+ \right] \right\}.
\end{eqnarray*}
Similarly, from the equality $H_1(x)=H_2^{\,'}(x)$ can be obtained $H_1(x)$. We compute that $x_1=-1-\sqrt{5}+2 \sqrt{2}$, $x_2=0$, $x_3=1+\sqrt{5}-2 \sqrt{2}$ are the points of sign changes of the function $H_2(x)$. It is not difficult to check that the assumptions of Corollary \ref{cor:2.7b} are satisfied.  Since
$$
(-1)^{3+1}H_3(x_2)=(-1)^{3+1}H_3(0)= \tfrac{1}{72}+\tfrac{\sqrt{5}}{360}-\tfrac{\sqrt{2}}{72}>0,
$$
it follows that the inequalities \eqref{eq:2.7d} are satisfied. From Corollary \ref{cor:2.7b} we conclude that the relation \eqref{eq:w4} hold.

\bibliographystyle{elsarticle-num}


\section{Introduction}

The approach to producing unconditional simulations for arbitrary surfaces (sphere, ellipsoid, geoid, torus, etc.) will be described in this paper. The approach is also applicable to nonsmooth surfaces. This paper is based on ideas from a turning band approach (\cite{TwoDimensionalSimulationByTurningBands}, \cite{TurningBandTwoDimensional}, \cite{ClosedFormSolutionsOfTheTwoDimensionalTurningBandEquation}, and \cite{TurningBandSimulationsInThreeDimensions}) and kernel convolution (\cite{KernelConvolution1996}, \cite{KernelConvolution1998}, and \cite{KernelConvolutionFFT}). Similar to the turning band, all processing is done on a set of random lines. The difference is in using a random set of points covering an area of interest and projecting them onto the set of random lines. For example, the set of uniformly distributed points on the surface of a sphere is projected onto the set of random lines. The kernel convolution is applied to each line, and the reconstruction of the process is performed similarly to the turning band approach. Using different kernels will produce a nonstationary solution.

\section
{
  Kernel Convolution by Turning Band
  \label{sec:KernelConvolutionByTurningBand}
}

Here are the steps of the algorithm:
\begin{enumerate}
  \item \label{en:UniformDistributionOnSurface} From a uniform distribution, simulate $N$ points $\vec{s_i}$ on a surface. \footnote{Replacing uniform distribution on the surface by pseudorandom uniform distribution will produce an even coverage and guarantee the stability of the algorithm.}
  \item From a standard normal distribution, simulate value $\xi_i$ for each point.
  \item Simulate $N_d$ random or pseudorandom directions (the algorithm to simulate pseudorandom directions can be found in \cite{RandomDirections} and \cite{TurningBandSimulationsInThreeDimensions}).
  \item Cycle over all directions.
  \begin{enumerate}
    \item Project points to a direction.
    \item Apply kernel convolution (this is usually done using FFT).
  \end{enumerate}
  \item For each location where the result is desired, project that location to all directions and average the results of kernel convolutions.
\end{enumerate}

The applied one-dimensional kernel $k_1{\left( \bullet \right)}$ will correspond to a three-dimensional kernel $k_3{\left( \bullet \right)}$ through integral
\begin{equation}
  k_3{\left( \vec{h} \right)}
  =
  \oiint
  {
    k_1{\left( \left( \vec{h}, \hat{n} \right) \right)}
    \dif{\hat{n}}
  }
  ,
  \label{eq:IntegralKernel}
\end{equation}
where $\left( \vec{x}, \vec{y} \right)$ is a scalar product.

The kernel distance is the direct distance in three-dimensional space.

This is similar to the turning band method, and the relationship between three-dimensional or two-dimensional covariance with one-dimensional covariance can be found in \cite{TwoDimensionalSimulationByTurningBands}, \cite{TurningBandTwoDimensional}, and \cite{ClosedFormSolutionsOfTheTwoDimensionalTurningBandEquation}. The same relationship holds true for the kernel functions. In a three-dimensional case, the formula is
\begin{equation}
  k_1(h) = (h \cdot k_3{\left( h \right)})'
  .
  \label{eq:RelationshipBetweenKernelsThree}
\end{equation}

In two dimensions for kernel $k_2{\left( \bullet \right)}$
\begin{equation*}
  k_1(h)
  =
  k_2{\left( 0 \right)}
  +
  h
  \int\limits_{0}^{h}
  {
    \dfrac
    {
      k_2'{\left( x \right)}
    }
    {
      \left(
        h^2 - x^2
      \right)
      ^
      {
        \frac{1}{2}
      }
    }
    \dif{x}
  }
  =
  k_2{\left( 0 \right)}
  +
  h
  \int\limits_{0}^{\frac{\pi}{2}}
  {
    k_2'
    {
      \left(
        h
        \sin{\left( \Theta \right)}
      \right)
    }
    \dif{\Theta}
  }
  .
\end{equation*}

\section{Kernel Convolution for Unconditional Simulation}

Applying the algorithm from the section ``\nameref{sec:KernelConvolutionByTurningBand}'' for any location $\vec{s}$ will produce a random field approximately equal to
\begin{equation}
  \sum_{i = 0}^{N}
  {
    \left(
      k_3
      {
        \left( \vec{s_i} - \vec{s} \right)
      }
      \cdot
      \xi_i
    \right)
  }
  .
  \label{eq:KernelApplied}
\end{equation}

The approximation is due to the finite number of directions.

The variance at location $\vec{s}$ is approximately equal to
\begin{equation}
  \sum_{i = 0}^{N}
  {
    k_3
    ^2
    {
      \left( \vec{s_i} - \vec{s} \right)
    }
  }
  .
  \label{eq:KernelVariance}
\end{equation}

For the sphere, the estimate of \eqref{eq:KernelVariance} does not depend on location $\vec{s}$. This does not hold true for other surfaces. However, for any surface with finite sampling $N$, it is not a constant. Calculating this value is important in standardizing unconditional simulations. The estimate of \eqref{eq:KernelVariance} can be found by applying the algorithm from the section ``\nameref{sec:KernelConvolutionByTurningBand}'' for the square of the kernel for the same set of points with values of one. Dividing the approximation of \eqref{eq:KernelApplied} by the square root of the approximation of \eqref{eq:KernelVariance} will produce an approximation of the Gaussian process with unit variance.

Another way to estimate \eqref{eq:KernelVariance} is by producing many unconditional simulations and estimating the variance from them. This will also allow the estimation of covariances.

Compared to the turning band approach, using this algorithm is significantly slower. The slowest part is projecting simulated points to directions. Notice that the simulated points do not depend on the kernel. Therefore, to speed up the approach, it is possible to simulate points and project them to directions in advance. Combining them by bins will significantly reduce the amount of memory necessary to store precalculated results.

\section{Conditional Simulations}

To perform conditional simulations, it is necessary to know the covariance function. For a sphere, the covariance function is isotropic and does not depend on location. Therefore, it can be derived directly from the kernel as described in \cite{CompactCovarianceModelOnASphere} and \cite{KernelConvolutionForRingsOnPlane}. For arbitrary surfaces, it is possible to make many unconditional simulations and from them estimate the covariance function. Conditional simulations are obtainable from unconditional simulations through solutions of the Kriging system \cite{ConditionalSimulationsGeostatisticsForConditionalSimulationOfOreBodies}, \cite{ConditionalSimulationsBook}, and \cite{ConditionalSimulationsSpatialPredictionSpatialSamplingAndMeasurementError}. \cite{GeostatisticalMappingWithContinuousMovingNeighborhood} describes the approach to construct conditional simulations without abrupt changes.

\section{Class of Exponential Kernels}

Let the kernel function have the form
\begin{equation}
  k_3{\left( h \right)}
  =
  \sum_{i}
  {
    \left(
      P_i{\left( h \right)}
      \cdot
      \mathsf{e}^{-\Theta_i \cdot h}
    \right)
  }
  ,
  \label{eq:Form}
\end{equation}
where $P_i$ is a polynomial.

The result of the convolution can be obtained by the sequential application of the kernel. The technique is described in \cite{Smoothing} and \cite{SmoothingNetwork} with the difference being that the integration is replaced by summation. This reduces the complexity of the kernel convolution from $ O{\left( n \cdot \ln{n} \right)} $ to $ O{\left( n \right)} $, where $n$ is the number of bins on the line.

For finite kernel
\begin{equation}
  k_3{\left( h \right)}
  =
  \left\{
  {
    \begin{aligned}
      &\sum_{i}
      {
        \left(
          P_i{\left( h \right)}
          \cdot
          \mathsf{e}^{-\Theta_i \cdot h}
        \right)
      }
      ,
      &\text{ if }
      h < T,\\
      &0,
      &\text{ otherwise},
    \end{aligned}
  }
  \right.
  \label{eq:FormFinite}
\end{equation}
where $T$ is the kernel width. To be able to evaluate $k_1{\left( h \right)}$ from \eqref{eq:RelationshipBetweenKernelsThree}, it is necessary to satisfy $\lim\limits_{h \rightarrow T^{-}}{k_3{\left( h \right)}} = 0$. Otherwise, $k_1{\left( h \right)}$ will not be finite.

There are several advantages of using the finite kernel:
\begin{itemize}
  \item Ability to produce unconditional simulations for a large area by partitioning it into small pieces (tiles) for independent processing.
  \item Ability to have a flexible class of kernels in form \eqref{eq:FormFinite}.
\end{itemize}

The ability to have a flexible class of kernels is important for nonstationary simulations. It allows not only changes in kernel bandwidth but in their shape as well. One way to construct a flexible class of kernels in form \eqref{eq:FormFinite} is by changing the polynomial part. It is reasonable to use decreasing kernels because the majority of physical processes gradually decreases dependence over distance. From~\cite{MonotonicFunction}, a Bernstein density (a weighted sum of beta probability distribution functions defined on the unit interval
$ \left[ 0 , 1 \right] $) is defined as
\newcommand{\myfunctionbeta}[1]
{
  Beta
  {
    \left(
      #1
      \mid
      j , m - j + 1
    \right)
  }
}
\newcommand{\myfunctionbetacorrection}[2][]
{
  \sum _{j = 1} ^{m}
    \left(
      \myfunctionwi[#1]{j}
      \cdot
      \myfunctionbeta{#2}
    \right)
}
\newcommand{\myfunctionbetacorrectionintegrated}[2][]
{
  \sum _{j = 1} ^{m}
    \left(
      \myfunctionwi[#1]{j}
      \int\limits_{x}^{1}
      {
        \myfunctionbeta{#2}
        \dif{y}
      }
    \right)
}
\begin{equation*}
  p{\left( x \mid m , \vec{w} \right)} = \myfunctionbetacorrection{x} \; , \;
  \myfunctionwi{j} \geq 0 \; , \;
  \sum _{j = 1} ^{m} {\myfunctionwi{j} = 1} \; ,
\end{equation*}
where
$ m $ is the number of distribution functions,
$ \vec{w} $ are weights,
$ Beta{\left( \bullet \right)}$
is the beta probability density function
\begin{equation*}
  \myfunctionbeta{x}
  =
  \dfrac{m !}{\left( j - 1 \right) ! \cdot \left( m - j \right) !}
  \cdot
  x ^{j - 1}
  \cdot
  \left( 1 - x \right) ^{m - j}
  \; .
\end{equation*}

Forming
\begin{equation*}
  P{\left( x \right)}
  =
  \int\limits_{x}^{1}
  {
    p{\left( x \mid m , \vec{w} \right)}
  }
  =
  \myfunctionbetacorrectionintegrated{y}
\end{equation*}
and multiplying by $\mathsf{e}^{-\frac{h}{\Theta}}$ produces a flexible family of kernels. The choice of $\Theta$ is arbitrary. If the value is too large, the floating point evaluation will lead to a high round off error. If the value is too small, it will not produce good family of kernels. Let's define $\Theta = 3$.

In the case of $m = 2$, this becomes
\begin{equation}
  \left\{
    \begin{aligned}
      &
      \left(
        \myfunctionwi{1} \left( 1 - h \right)^2
        +
        \myfunctionwi{2} \left( 1 - h^2 \right)
      \right)
      \cdot
      \mathsf{e}^{-\frac{h}{3}}
      ,
      &
      \text{ if }
      h < 1,
      \\
      &
      0
      ,
      &
      \text{otherwise},
    \end{aligned}
  \right.
  \label{eq:KernelsPolynomialExponentTwo}
\end{equation}
where $0 \leq \myfunctionwi{i} \leq 1$, $i = 1, 2$, and $\myfunctionwi{1} + \myfunctionwi{2} = 1$. They are shown in Figure~\ref{fig:KernelsPolynomialExponent}a.

\begin{figure} [htb]
  \centering
  \begin{tabular}{c c}
    \includegraphics[width = 8 cm, keepaspectratio]{KernelsPolynomialExplonent.png} &
    \includegraphics[width = 8 cm, keepaspectratio]{KernelsPolynomial.png} \\
    a) &
    b)
  \end{tabular}
  \caption
  {
    Example of kernels \eqref{eq:KernelsPolynomialExponentTwo}.
    $\myfunctionwi{1} = \dfrac{i}{15}$, $\myfunctionwi{2} = 1 - \myfunctionwi{1}$, $i = \overline{0..15}$. a) With exponent. b) Without exponent.
  }
  \label{fig:KernelsPolynomialExponent}
\end{figure}

In the case of $m = 3$
\begin{equation}
  \left\{
    \begin{aligned}
      &
      \left(
        \myfunctionwi{1} \left( 1 - h \right)^3
        +
        \myfunctionwi{2} \left( 1 - h \right)^2 \left( 1 + 2 h \right)
        +
        \myfunctionwi{3} \left( 1 - h^3 \right)
      \right)
      \cdot
      \mathsf{e}^{-\frac{h}{3}}
      ,
      &
      \text{ if }
      h < 1,
      \\
      &
      0
      ,
      &
      \text{otherwise},
    \end{aligned}
  \right.
  \label{eq:KernelsPolynomialExponentThree}
\end{equation}
where $0 \leq \myfunctionwi{i} \leq 1$, $i = 1, 2, 3$, and $\myfunctionwi{1} + \myfunctionwi{2} + \myfunctionwi{3} = 1$. They are shown in Figure~\ref{fig:KernelsPolynomialExponentThree}a.

\begin{figure} [htb]
  \centering
  \begin{tabular}{c c}
    \includegraphics[width = 8 cm, keepaspectratio]{KernelsPolynomialExplonentThree.png} &
    \includegraphics[width = 8 cm, keepaspectratio]{KernelsPolynomialThree.png} \\
    a) &
    b)
  \end{tabular}
  \caption
  {
    Example of kernels \eqref{eq:KernelsPolynomialExponentThree}.
    $\myfunctionwi{1} = \dfrac{i}{7}$, $\myfunctionwi{2} = \dfrac{j}{7}$, $\myfunctionwi{3} = 1 - \myfunctionwi{1} - \myfunctionwi{2}$, $i = \overline{0..7}$, $j = \overline{0..7 - i}$. a) With exponent. b) Without exponent.
  }
  \label{fig:KernelsPolynomialExponentThree}
\end{figure}

From \eqref{eq:RelationshipBetweenKernelsThree}, it follows that $k_1{\left( h \right)}$ has the same form as \eqref{eq:FormFinite}. Notice that $k_3^2{\left( h \right)}$ also has the same form, and the one-dimensional kernel corresponding to $k_3^2{\left( h \right)}$ will be of the same form.

It is possible to extend this approach for anisotropic kernels by adjusting weights and one-dimensional kernel shape in each direction or by stretching the space.

\section{Examples}

Figure~\ref{fig:UnconditionalSimulationTorus} shows unconditional simulation on a torus.

\begin{figure} [htb]
  \centering
  \includegraphics[width = 12 cm, keepaspectratio]{UnconditionalSimulationTorus.jpg}
  \caption{Unconditional simulation on a torus with inner radius 1, outer radius 3, and kernel $\mathsf{e}^{-10 h}$. Number of simulated points $N = 67,108,864$. Number of random directions $N_d = 1,024$.}
  \label{fig:UnconditionalSimulationTorus}
\end{figure}

Figure~\ref{fig:UnconditionalSimulatiOnEllipsoid} shows nonstationary unconditional simulation on a spheroid with eccentricity equals $\frac{3}{4}$. The kernel changes smoothly in shape from left to right. $1,024$ unconditional simulations were used to estimate variance and standardize unconditional simulations. Due to changes in the kernel, the right and left parts of the spheroid have different covariance structures. This is visible as a smoother random field on the right compared to the random field on the left.

\begin{figure} [htb]
  \centering
  \includegraphics[width = \columnwidth, keepaspectratio]{NonstationarySimulationsSpheroid.png}
  \caption
  {
    Unconditional simulation on an oblate spheroid with semimajor axis equals $2$, and semiminor axis equals $1$. The surface is shown as elevation to depict changes in smoothness. The kernel has the form of \eqref{eq:KernelsPolynomialExponentTwo}. It changes smoothly from the kernel on the left
    $
      \left( 1 - \left( 3h \right)^2 \right)
      \cdot
      \mathsf{e}^{-h}
    $
    to the kernel on the right
    $
      \left( 1 - 3h \right)^2
      \cdot
      \mathsf{e}^{-h}
    $
    with bandwidth equals $\frac{1}{3}$.
    Number of simulated points $N = 262,144$. Number of random directions $N_d = 1,024$.
  }
  \label{fig:UnconditionalSimulatiOnEllipsoid}
\end{figure}

\section{Convergence Issue}

When the kernel bandwidth tends to be small, the convergence of \eqref{eq:KernelApplied} to the true kernel \eqref{eq:IntegralKernel} becomes slow. The number of required directions makes this method inefficient. This is because the distribution of the number of points that are farther away from the center of the kernel (location $\vec{s}$ in \eqref{eq:KernelApplied}) in three dimensions increases squarely. This is even true for some directions when points are only located on the surface.

\section{Tile Approach}

Finite kernel and overlapping tiles overcome the issue described in the previous section. For example, in two dimensions, if all tiles are $2$ by $2$ and start from integer numbers (corners can be rounded), then the kernel with a bandwidth not exceeding $\frac{1}{2}$ will always be completely inside at least one tile; see Figure~\ref{fig:OverlappingTiles}. This is applicable in any number of dimensions.

\begin{figure} [htb]
  \centering
  \begin{tikzpicture}[scale = 2.0]

    \fill [cyan] (2.5, 1) -- (3.5, 1) arc (-90:0:0.5) -- (4, 2.5) arc (0:90:0.5) -- (2.5, 3) arc (90:180:0.5) -- (2, 1.5) arc (180:270:0.5) -- cycle;

    \fill [red] (2.7, 1.8) circle (0.5);
    \tkzInit[xmin = 0, ymin = 0, xmax = 5, ymax = 4]
    \tkzGrid
    \tkzAxeXY
  \end{tikzpicture}
  \caption
  {
    The red circle is the area covered by a finite kernel. The blue rounded square is the tile completely covering the finite kernel.
  }
  \label{fig:OverlappingTiles}
\end{figure}

The bandwidth should not be significantly smaller than the tile width, and the exponential part in the kernel can be dropped. The flexible family of kernels without the exponential part is shown in Figures~\ref{fig:KernelsPolynomialExponent}b and \ref{fig:KernelsPolynomialExponentThree}b.

Figure~\ref{fig:UnconditionalSimulatiOnPlane}a shows the effect of using the tile approach. The small discontinuities are visible in the horizontal and vertical lines passing through the center. Increasing the number of bands will reduce this effect. Figure~\ref{fig:UnconditionalSimulatiOnPlane}b shows the result of evaluating \eqref{eq:KernelApplied} and \eqref{eq:KernelVariance} directly.

\begin{figure} [htb]
  \centering
  \begin{tabular}{c c}
    \includegraphics[width = 8.5 cm, keepaspectratio]{UnconditionalSimulationsUsingTiles.png} &
    \includegraphics[width = 8.5 cm, keepaspectratio]{UnconditionalSimulationsExact.png} \\
    a) &
    b)
  \end{tabular}
  \caption
  {
    Unconditional simulation in two dimensions by kernel
    $
      k_3{\left( h \right)}
      =
      \left( 1 - 2h \right)^2
    $
    with bandwidth equals $\frac{1}{2}$.
    Number of simulated points $N = 65,536$. a)~The tile approach with the number of random directions $N_d = 1,024$ (uniformly distributed in three dimensions). b)~Direct calculation.
  }
  \label{fig:UnconditionalSimulatiOnPlane}
\end{figure}

If a continuous result is required, overlapped tiles can be joined smoothly. This will increase calculations by the number of tiles necessary to process. One way to reduce the number of tiles is to use a triangular grid (use nonequilateral triangles for arbitrary surfaces). For each node of the triangular grid, define a tile that covers all neighboring triangles and at least a kernel bandwidth. For each triangle, the barycentric coordinate system will give weights for the tiles. There will only be a maximum of three tiles with nonzero weights; therefore, the complexity is increased three times.

\section{Using Integer Arithmetic}

In the turning band approach, the directions are chosen randomly or pseudorandomly. In this section, another approach will be described to make directions proportional to integer numbers. Using such directions will remove errors from applying fixed steps (bins) in each direction.

Let's take a sphere of radius $R = \sqrt{S}$, where $S$ is a natural number. If the radius $R$ approaches infinity, all integer vectors inside the sphere (all integer vertices except the center of a coordinate system) tend to be uniformly distributed. However, the unique directions are not uniformly distributed (see example in Figure~\ref{fig:IntegerDirections}a). But assigning a proper weight to each direction will resolve it. The appropriate weights can be assigned as proportionate to the number of vectors with the same direction or to the area of the spherical Voronoi diagram cells.

\begin{figure} [htb]
  \centering
  \begin{tabular}{c c}
    \includegraphics[width = 8.5 cm, keepaspectratio]{IntegerDirections.png} &
    \includegraphics[width = 8.5 cm, keepaspectratio]{IntegerUniformDirections.png} \\
    a) &
    b)
  \end{tabular}
  \caption
  {
    a) All integer directions are shown as dots on the sphere for $S = 100$. The number of unique directions equals $1,729$ (counting $\vec{v}$ and $-\vec{v}$ once). b) A uniform set of integer directions satisfying that any two directions are at least $3 \degree$ apart. The number of directions equals $1,405$.
  }
  \label{fig:IntegerDirections}
\end{figure}

Selecting only directions separated by some degree will make a uniform coverage (see Figure~\ref{fig:IntegerDirections}b). The list of integer directions is shown in the \hyperref[SectionAppendix]{appendix}. This is similar to using a set of pseudorandom directions as described in \cite{RandomDirections} and \cite{TurningBandSimulationsInThreeDimensions}.

Because each direction is represented as $\alpha \cdot \left( x, y, z \right)$, where $\alpha$ is a real number (square root of a natural number) and $x$, $y$, and $z$ are integers, any coordinate with integer components will be projected to a position $\alpha \cdot i$, where $i$ is an integer number.

If random points are placed on an integer grid, there are no errors due to binning. If the resultant location also has an integer coordinate, there is no additional error. Notice that the kernel \eqref{eq:Form} or \eqref{eq:FormFinite} can be calculated for any position on direction, see \cite{Smoothing} and \cite{SmoothingNetwork}. The turning band approach \cite{TwoDimensionalSimulationByTurningBands}, \cite{TurningBandTwoDimensional}, \cite{ClosedFormSolutionsOfTheTwoDimensionalTurningBandEquation}, and \cite{TurningBandSimulationsInThreeDimensions} might also benefit from using integer arithmetic.

\section{Conclusion}

The new approach to producing unconditional nonstationary simulations on arbitrary surfaces is equivalent to the kernel convolution approach with the kernel shape depending on location (\cite{KernelConvolution1996}, \cite{KernelConvolution1998}, and \cite{KernelConvolutionFFT}). Applying a kernel to the arbitrary surfaces produces a larger class of valid covariance functions than the class of valid covariance functions in three-dimensional space.

Precalculating projections of random points in all directions for all tiles might be necessary for the efficiency of the approach. Nevertheless, even with precalculations, this approach is not efficient. The better approach is based on kernel convolution using FFT as described in \cite{ConvolutionFFT}.

New flexible classes of kernels constructed from polynomials are described in this paper. Because they are linear combinations of some basis kernels, constructing covariances for all pairs of basis kernels is sufficient to calculate covariance between any two kernels. Unconditional simulation for a nonstationary kernel is a weighted sum of unconditional simulations corresponding to basis kernels for the same set of random points.

Using integer arithmetic, an error introduced in the turning band approach due to binning along each direction can be avoided.

\section{Acknowledgment}

The author would like to thank David Burrows for help with area preserving transformation from authalic sphere to spheroid (authalic sphere is a sphere were the surface area equals the spheroid surface area). This transformation was used in simulating uniformly distributed random points on the surface of a spheroid for the creation of Figure~\ref{fig:UnconditionalSimulatiOnEllipsoid}.

\newcommand{\doi}[1]{\textsc{doi}: \href{http://dx.doi.org/#1}{\nolinkurl{#1}}}

\begingroup
\raggedright
\bibliographystyle{IEEEtran}


\section{Introduction}
The Southern hemisphere site of the Cherenkov Telescope Array (CTA) \cite{CTA} is envisaged to host $\sim$70 Small-Sized Telescopes (SST) \cite{SST}, optimised for the highest gamma-ray energies that CTA will be sensitive to. The Gamma-ray Cherenkov Telescope (GCT), a sub-consortium of CTA, proposes to construct 35 SSTs, designed to cover the energy range from about 1 to 300 TeV. The GCT telescope is a dual mirror Schwarzschild-Couder design with a primary mirror of diameter 4~m, a secondary mirror of diameter 2~m, a focal length of 2.3~m and a focal ratio of 0.57 (see \cite{GCT_Telescope}) resulting in a spherical focal plane with a radius of curvature 1.0~m. The GCT camera is designed to record flashes of Cherenkov light lasting from a few to a few tens of nanoseconds, with typical image width and length of $\sim$0.2$^{\circ} \times $1.0$^{\circ}$ and promises a low-cost, high-reliability, high-data-quality solution for a dual-mirror SST. The small focal length of the telescope implies that an approximately 0.2$^{\circ}$ angular pixel size (as required by CTA for the SSTs) is achievable with pixels of physical dimensions of 6 to 7 mm, while the dual-mirror optics ensure that the point spread function (PSF) of the telescope is below 6~mm up to field angles of 4.5$^{\circ}$.  A field of view (FoV) of 8$^{\circ}$ (again required by CTA for the SSTs) can therefore be covered with a camera of diameter about 0.4 m, composed of 2048 pixels. This allows the use of commercially available photosensor arrays, significantly reducing the complexity and cost of the camera. Suitable packages consisting of multi-anode photomultipliers (MAPMs) or silicon photomultipliers (SiPMs) are under investigation for the GCT via the development of two prototypes (CHEC-M) and (CHEC-S) respectively, where CHEC stands for Compact High Energy Camera. In this paper we focus on the former which is now under testing in the lab. In section \ref{GCT-SiPM} we briefly report the ongoing work towards the completion of the second camera.

\section{CHEC-M Design}

A CAD image of CHEC-M is shown in Figure~\ref{FigCam} showing the principle external components, whilst Figure \ref{FigElectronics} shows the internal electronics architecture for the camera. The camera is designed to provide full waveform information for every pixel in every event. It is based on 32 camera modules each containing a 64-channel MAPM attached via a preamplifier board to a front-end module called TARGET providing digitisation and triggering. The modules are mounted in an internal rack and directly to a backplane PCB, which in-turn routes raw data to two data acquisition boards (DACQ). The position of the MAPMs at the focal plane is defined and aligned on the focal plane with the required accuracy of better than 0.35 ~mm. Individually shielded ribbon cables between the preamplifiers and TARGET modules are used to match the different physical geometries of the internal rack and the curved focal plane. LED flasher units placed in the corner of the camera provide calibration over a range of illumination intensities via reflection from the secondary mirror (see \cite{Flasher} for further details). Internal fans circulate air via a set of mechanical baffles. An external lid system provides protection from the elements. The lid, fans, LED flashers and a set of internal sensors (temperature and humidity) are controlled via a peripherals board. The camera is liquid cooled via a chiller unit mounted on the telescope. Camera readout and communications take place via optical fibre. The camera is powered by a single 12~V DC supply, mounted at the rear of the secondary mirror, and the power consumption is $\sim$450~W.

\begin{figure}[!ht]
  \centering
\resizebox{1\columnwidth}{!}{\includegraphics[trim=1cm 0cm 1cm 0cm, clip=true]{GCT_Cam.pdf}}
\caption{A detailed overview of the GCT camera. The CAD model shown is
  for the prototype camera equipped with MAPMs (CHEC-M). In addition a protective window
  may be used to protect the photodetectors.}
  \label{FigCam}
\end{figure}

\begin{figure}[!ht]
  \centering
  \makebox[1\columnwidth][c]{\includegraphics[trim=3.5cm 1.5cm 4.5cm 0cm,
    clip=true, width=1.16\textwidth]{GCT-Cam-Architecture_v2.pdf}}
  \caption{A schematic showing the logical elements of the GCT camera,
    the communication between those elements, the raw data flow
    through the camera, the trigger architecture and the clock
    distribution scheme. Power distribution is excluded for
    simplicity.}
  \label{FigElectronics}
\end{figure}

\subsection{Front-End Electronics}

Figure~\ref{FigTM} shows the front-end camera module used in CHEC-M. Each MAPM has 64 pixels with physical size of $\sim$6~mm$\times$6~mm. The photosensor is connected to a preamplifier module, containing 4 PCBs of 16 channels each, which shape the signal from the fast detector stretching the signal to allow for an optimal coincidence time window for the trigger. This optimal shape is determined from Monte Carlo simulations to be: 5.5--10.5~ns FWHM, 3.5--6.0~ns rise time. The shaped pulses are then fed to the TARGET module via ribbon cables. The TARGET module is based on four custom TARGET ASICs \cite{TARGET}, each one of which accepts 16 analogue inputs from the preamplifier module. The shaped signal is digitised with 12-bit resolution at 1~GS/s over a programmable read out window in 32~ns steps (96~ns is currently envisaged for normal operation). The conversion and readout lead to a dead time of $\sim$20 $\upmu$s, which is acceptable at the maximum average camera trigger rate of $\sim$600~Hz. The CHEC-M TARGET modules use TARGET 5 ASICs which show an RMS noise of 0.5~mV and a dynamic range of 1.2~V. Resulting in a pixel saturation at $\sim$500 p.e. assuming operation at a gain level such as 1~p.e. corresponds to 5 times the noise RMS. Simulations have shown that the waveform information can be used to extend the dynamic range beyond saturation. The TARGET module also provides the first level of trigger based on the discrimination of the analog sum of 4 neighbouring pixels (super-pixel). The TARGET module contains a Xilinx Spartan-6 FPGA to process trigger signals (dark and light blue lines in Fig. \ref{FigElectronics}) and read waveforms(gray line in Fig. \ref{FigElectronics}) from arbitrary storage capacitors. Furthermore it provide interface for slow control: set the supplied HV to the MAPM (operating range 800--1100 V), monitor status of the module, configure and read out the ASICs. Serial data is then output from the TARGET module over ethernet using UDP protocol for raw data, control and monitoring.

\begin{figure}[!ht]
\begin{centering}
\resizebox{1\columnwidth}{!}{\includegraphics[trim=3cm 3cm 3cm 2cm, clip=true]{TM.pdf}}
\caption{Photograph of a TARGET 5 based front-end camera module used in CHEC-M.}
\label{FigTM}
\end{centering}
\end{figure}

\subsection{Back-End Electronics}
A single backplane PCB is used to gather the signals from all the 32 TARGET modules. The backplane routes the UDP raw data lines to the DACQ boards, distributes a synchronous clock signal, supplies the TARGET modules with power (and control of that power), and provides the camera-level trigger. A Xilinx Virtex 6 FPGA accepts all 512 LVDS first-level trigger lines from the TARGET modules and implements a camera-level trigger algorithm (currently any two triggered neighbouring super pixel patches trigger the camera). Upon a valid trigger the TARGET modules are read out and a unique event number is added to the raw data. In the final CTA-scheme a timing board based on the White Rabbit protocol\cite{WhiteRabbit} will be used to provide an absolute timestamp and associate it with the unique event number. For CHEC-M the DACQ boards are capable of implementing a White Rabbit timing interface on-board. All raw data, control and monitoring communication to the camera is via four 1~Gbps fibre optic cables attached to the DACQ boards.

\section{CHEC-M Commissioning}
Figure~\ref{FigChecM} shows the fully assembled prototype camera populated with all front-end modules, internal mechanics and DACQ boards. A prototype backplane PCB is used to verify camera readout whilst the real backplane is undergoing separate lab-tests and firmware development. The camera is held horizontally on a sliding rail in a dark box also hosting a robot arm holding a near UV laser as shown in Figure~\ref{FigChecM}. For the purpose of flat fielding a diffuser is mounted in front of the laser and the robot arm is held as far from the camera as the physical size of the dark box allows to have the most uniform light beam possible.

Calibration and refining of the data-taking process is in progress. Figure~\ref{FigCom1} shows some initial commissioning results from CHEC-M, showing the pulse area as measured from an approximately uniform light flash (top, middle) and example waveforms from the same data for an individual pixel (top, right, bottom) and averaged over an MAPM (top, right, top). The bottom images are taken with a mask in front of the camera to approximate a Cherenkov image with a cut on the pulse peak height. In both cases, minimal gain-matching between MAPMs has been applied and pedestals have been subtracted, but no other calibration has been applied.

Once commissioned, testing will include illuminating the full camera with uniform flashes of light at varying brightness (using both a laser and a CHEC LED flasher) in the presence of a realistic background light level (from a diffuse white-light source). The robot arm will also allow the focal plane to be scanned with a narrow beam from the laser to measure the variation in response as the edge of pixels and gaps between MAPMs are encountered as well as the angular response of the pixels.

\begin{figure}
\begin{centering}
\resizebox{1\columnwidth}{!}{\includegraphics[trim=2cm 6.5cm 3cm 0cm, clip=true]{CHEC-M.jpg}}
\caption{CHEC-M: (a) MAPM photodetectors, (b) preamplifier modules, (c)
  TARGET modules, (d) proto-Backplane, (e) FPGA test board, (f) DACQ
  boards, (g) calibration flashers, (h) internal mechanics consisting
  of focal plane plate and electronics rack, (i) 12 V input cables (j)
  lab fan-cooling, (k) sliding rail, (l) illumination source, (m) robot arm, (n) dark box, (o) external mechanical enclosure.}
\label{FigChecM}
\end{centering}
\end{figure}

The GCT camera has to operate reliable in a very harsh environment for more than 15 years with only minor maintenance. Therefore an empty enclosure was artificially exposed in the lab to verify performance against CTA requirement. The tests performed were: impact testing, water ingress, temperature cycling (followed by further water-ingress testing), UV exposure and wind testing. Wind-tunnel tests showed we can operate the lid up to a maximum wind load of 50~km/hour.  Impacts on both the camera body and lid showed no ill effects, the camera body (including the interface with the faceplate) showed no evidence of water ingress, and there was no damage during temperature cycling between ${-}25$ and +40$^\circ$C or during UV exposure. Water ingress was found to be present via the seal between the camera lid and body. This lead to a minor, but important redesign to ensure camera survival on site.

On-telescope tests are planned to take place in Paris on the GCT Telescope prototype structure. On-telescope testing will include: the alignment of the camera in the focal plane, in-situ calibration, recording Cherenkov images, and procedural aspects of installing, operating and maintaining the camera as part of the telescope system.

\begin{figure}
\begin{centering}
\resizebox{0.9\columnwidth}{!}{\includegraphics[trim=2.2cm 7cm 2.5cm 1cm,
  clip=true]{CHEC-M_Commissioning1.pdf}}
\resizebox{1.0\columnwidth}{!}{\includegraphics[trim=1cm 11cm 2cm 1cm, clip=true]{CHEC-M_Commissioning2.pdf}}
\caption{Initial commissioning results from CHEC-M measuring uniform light flashes. (top, middle) pulse
  area, (top, right, top) example waveforms for an individual
  pixel and (top, right, bottom)  averaged of an MAPM. (bottom) images taken with a mask in-front of the camera
  to approximate the shape of Chernekov light flashes from gamma-ray initiated showers.}
\label{FigCom1}
\end{centering}
\end{figure}

\section{Beyond CHEC-M}\label{GCT-SiPM}
The second prototype GCT camera (CHEC-S) features similar mechanics and the same back-end electronics as CHEC-M. SiPMs are used instead of MAPMs and an upgrade of ASIC from TARGET 5 to TARGET 7 is present. A revised pre-amplification and shaping scheme is needed to accommodate the SiPM pulse shape and an update to the TARGET modules is implemented. The TARGET 7 ASICs feature improved linearity and dynamic range when compared to TARGET 5. CHEC-S is expected to be fully integrated towards the end of 2015. Beyond CHEC-S the GCT camera will go through one further design iteration before entering the CTA Production phase. During this design iteration the final photosensor choice will be made and a final iteration of the TARGET ASICs will take place.

\section{Summary}

The GCT camera is a low-cost, high-reliability, high-data-quality solution for the SST section of CTA. The camera features full waveform readout with nanosecond sampling across a variable window for all 2048 pixels at a rate of up to 600~Hz. The use of high-channel density TARGET ASICs and the latest advances in multi-pixel photosensor tiles ensure a cost-effective design. The first camera prototype (CHEC-M) has been assembled and is undergoing commissioning in the lab. After full characterisation the camera is scheduled to be installed on a prototype telescope structure in fall 2015. Beyond CHEC-M a second prototype based on SiPMs is under development with the eventual aim of providing 35 cameras to CTA.

\acknowledgments
The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Unions Seventh
Framework Programme FP7/2007-2013/ under REA grant agreement n [317446] INFIERI ``INtelligent Fast Interconnected and Efficient Devices
for Frontier Exploitation in Research and Industry``.We would also like to acknowledge the support of the UK Science and Technology Facilities Council (grant ST/K501979/1) and the support from the agencies and organisations listed in this page: http://www.cta-observatory.org.