\section*{Results}
\emph{Phase diagram --}
Using an approach similar to that used for solving the SK model, the exact $d\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle\io$ solution for
$d$-dimensional identical hard spheres of unit diameter
can be formulated in terms of a caging order parameter $\D(y)$~\cite{Pa83,MPV87}.
This functional order parameter, which encodes the width $\D$ of metabasins on a (properly defined) scale $y$,
is obtained by numerically solving a set of integro-differential equations (Supplementary Note~1)
and then used to calculate the theoretical liquid-glass phase diagram (Fig.~\ref{fig:PD}).
The theory predicts that a compressed liquid falls out of equilibrium and becomes a glass at a pressure that
depends on the compression rate. Once in a glass state, further compression results in a quick increase of the system pressure $p$, and upon jamming $p\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle \io$~\cite{PZ10}. The final jamming density depends
on compression speed, hence defining a J-line of jammed states~\cite{PZ10}. Two examples of glass compression obtained using an approximate state following are reported in Fig.~\ref{fig:PD}~\cite{KZ10} .

Independently of compression rate, the glass basin in which the system is initially trapped
undergoes a Gardner transition~\cite{Ga85}, at a line computed in Ref.~\cite{KPUZ13}.
Our key result is that, at pressures above this line, \emph{basins transform
into metabasins that contain a collection of marginally stable glasses},
a phenomenon that is described by a non-trivial caging order parameter $\D(y)$ as in the SK model~\cite{Pa83}.
Finding the solution that describes the marginal phase allows us to delimit the marginal phase boundary (Supplementary Note~1) to within the Gardner transition
line of Ref.~\cite{KPUZ13}, the J-line, and the ``threshold'' line determined following the prescription of Ref.~\cite{Ri13}.
The fact that within this region at least one eigenvalue of the stability matrix in the free energy space vanishes
confirms that this phase is indeed marginally stable (Supplementary Note~1)~\cite{DK83}. We also find that,
while the radius of the innermost fractal basins shrinks to zero as a power-law $\D_{\mathrm{EA}} \sim p^{-\k}$ (see below),
the radius of the largest metabasins remains of order one.
Close to jamming the total entropy of a group of metabasins of width $\Delta$ grows as $\Delta^{1/\k}$,
 hence the basins have a phase space structure whose fractal dimension is $2/\k$
(see Supplementary Note~1 for a more detailed discussion).
The {\it marginal} phase is thus also {\it fractal}.

The existence of the marginal phase can be qualitatively tested by molecular dynamics (MD) numerical simulations
in finite $d$ (Supplementary Note~2) by considering the outcome of a slow compression from the liquid up to
jamming~\cite{CIPZ11,CCPZ12}.
Jammed systems are isostatic,
and thus particles have an average of 2$d$ force-bearing neighbors~\cite{Mo98,TW99,Ro00,OLLN02,BW09b,CCPZ12},
which is much smaller than the $\mathcal{O}(e^d)$ neighbors that isotropically cage a particle in an equilibrated dense liquid.
Because the identity of the force-bearing neighbors at jamming uniquely characterizes the state,
their emergence sensitively depends on the landscape structure (Fig.~\ref{fig:basins}).
In the simple basin scenario, force-bearing neighbors at jamming should be fully determined immediately upon leaving the equilibrium liquid; in the meta/subbasin scenario, that determination should only occur once sufficiently deep in the glass for transitions between subbasins to be fully suppressed; in a fractal phase,
by contrast, the contacts should be gradually determined as jamming is approached.
To test this scenario, we consider a glass configuration at pressure $p_{\mathrm{init}}$. Starting from this configuration, we perform several
independent compressions up to $p_\mathrm{f} = 10^{10}$ and for each compressed configuration we measure the force network.
We obtain a set of contact variables $f^{(a)}_{ij}$, which are set to unity if particles $i$ and $j$
form a force-bearing contact in configuration $a$ and
to zero otherwise. The average of $\langle f^{(a)}_{ij} f^{(b)}_{ij} \rangle$ over pairs $ab$ of compressed configurations and over contacts $ij$ provides
a measure of similarity between the force networks. The fact that this quantity
increases smoothly upon increasing $p_\mathrm{init}$ indicates that
the force network is only partially encoded in the initial configuration, in support
of the fractal landscape scenario (Fig.~\ref{fig:MDSA}a and Supplementary Note~2).

\begin{figure}
\center{
\includegraphics[width=0.5\columnwidth]{MSDA.eps}}
\caption{
(a) Overlap between the force network edges $f_{ij}={0,1}$
that connect two particles $i$ and $j$ in two glass configurations $a$ and $b$ (at pressure $p_\mathrm{f} = 10^{10}$),
obtained by independent compression of the same initial configuration at initial pressure $p_\mathrm{init}$ (Supplementary Note~2).
(b) Time-evolution of the mean-square displacement $\Delta(t)$ for glasses at $p=10^{2}$, $10^3$, $10^4$, $10^5$, $10^6$, $10^7$, and $10^{8}$ in $d$=4. The solid line indicates the ballistic $dt^2$ behavior. The long-time value is the cage size $d \, \Delta_{\mathrm{EA}}$.
(c) The pressure evolution of the cage size $\Delta_{\mathrm{EA}}$ in  various dimension closely follows a power-law $\sim p^{-\k}$ with
$\k=1.41575$ as predicted by the theory.
(d) Analytical results for the order parameter $\D(y)$ at $\widehat} \newcommand{\wt}{\widetilde\f=10$.
Increasing the cutoff $y_{\rm max}\sim p$ indicates that the scaling regime $\D(y) \sim y^{-\k}$ extends to all
$y$. (Supplementary Note~1)}
\label{fig:MDSA}
\end{figure}

\emph{Criticality of the jamming transition --}
The equations that describe the marginal phase are formulated in terms of the caging order parameter $\Delta(y)$ and the pair correlation function $g(r)$, which also encodes the probability distribution of forces in the packing.
Upon approaching the J-line, i.e., as $p\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle\io$, these equations develop a scaling regime (Fig.~\ref{fig:MDSA}d) that is characterized by three main critical exponents: $\theta$ for the weak forces, $\alpha$ for the quasi-contacts, and $\kappa$ for $\D$ itself.
A (non-trivial) generalization of the approach developed for the SK model~\cite{Pa06} allows us to obtain theoretical values for these exponents (Supplementary Note~1). Interestingly, the condition that fixes their value is precisely equivalent to the marginal stability condition.
The theory therefore predicts that the criticality of the jamming transition directly follows from its location inside the marginal phase.

A striking signature of marginality is the scaling of the innermost basin width captured by the Edwards--Anderson cage size
$\D_{\mathrm{EA}} \sim p^{-\kappa}$. Although $\kappa=3/2$ was proposed in earlier studies~\cite{WNW05,WSNW05,BW09b,IBB13},
the theory predicts a slightly smaller $\kappa = 1.41574$ that is in remarkable agreement with our numerical results
(Fig.~\ref{fig:MDSA}c).
Because single-particle caging by immediate neighbors (a simple Einstein model for glasses),
would give $\k=2$~\cite{IBB13}, $\k < 2$ implies that fluctuations near jamming are divergently larger than for independent vibrations, in support of
their cooperative nature~\cite{BW09b,IBB13}.
Note that if one ignores the fractal phase, an explicit computation erroneously gives $\k=1$~\cite{KPUZ13}.
Also, note that the exponent $\k$ controls the fractal dimension of the basins, as discussed above.

The pair correlation function $g(r)$ also bears a signature of the criticality at the jamming transition.
The theory predicts, consistently with the analysis of~\cite{DTS05}, that when $p \rightarrow} \def\la{\left\langle} \def\ra{\right\rangle \io$, $g(r)$ develops an isostatic contact peak
characterized by a scaling function $\FF(\l) \equiv g(r)/g(1)$ for $\l=(r-1)p$ (Supplementary Note~1 and Fig.~\ref{fig:gr}).
It also predicts that the scaling
function of the contact peak decays as $\FF(\l)\sim \l^{-2-\th}$ at large $\l$.
The distribution $P(f)$ of inter-particle forces in the packing, which is related to the scaling function of the contact peak by
$\FF(\l)=\int_0^\io df \, f\, P(f) \, e^{-\l f}$~\cite{DTS05,PZ10,CCPZ12}, thus also decays as a power law $P(f)\sim f^\theta$
at small forces. Note that, as observed in~\cite{Wy12}, this phenomenon is closely related to what happens to the distribution of
frozen fields in the SK model~\cite{SD84},
which is thought to explain the Coulomb gap in interacting electron systems~\cite{MP07}.
Beyond the contact peak, the slower decay of pair correlation function follows another power-law $g(r) \sim (r-1)^{-\a}$ that describes the abundance
of quasi-contacts.
These scalings of $g(r)$ are crucial for determining the mechanical stability of packings~\cite{Wy12,LDW13}. Perturbing a packing breaks
some contacts with small forces, while also forming new contacts from what were previously quasi-contacts. Based on this observation a scaling relation for mechanical stability $\a = 1/(2+\th)$ can be derived~\cite{Wy12}. Remarkably, the exponents {\it predicted} by our theory, $\a=0.41269$ and $\th=0.42311$, satisfy this scaling relation to within numerical precision.
Prior estimates of these exponents were also obtained by numerical simulation. The quasi-contact exponent $\a$ has been measured by several groups in dimension $d$ ranging
from $2$ to $13$, all obtaining
roughly $\a\approx 0.4$~\cite{DTS05,SDST06,CCPZ12,LDW13,AST13}, the most precise estimates being $\a=0.41(3)$~\cite{CCPZ12}.
The weak force exponent $\th$ is, however, more difficult to measure,
and values spanning the interval $\th \approx 0.2 \div 0.45$ have been reported~\cite{CCPZ12,LDW13}.
Although the existence of a second exponent $\theta'<\theta$ has been shown to affect the tail of $P(f)$~\cite{LDW13},
its role in determining $\FF(\l)$ and its large-dimensional scaling remains to be clarified.
Additional numerical simulations are thus needed to test the theory more stringently.

A prediction for the force distribution $P(f)$ at jamming is also available from the theory.
But because this function is not completely determined by the scaling regime, it must be obtained by solving
the full equations that describe the marginal phase.
Numerically, the function $\FF(\l)=\int_0^\io df \, f\, P(f) \, e^{-\l f}$ is much easier to measure than $P(f)$ because it only depends on structural information, while, in
hard spheres, forces must be determined from the collision dynamics. A measure of $\FF(\l)$ thus also provides a more precise way to measure $\th$.
The theoretical prediction for the scaling function $\FF(\l)$ and therefore for $\th$ are
tested against numerical simulations in Fig.~\ref{fig:gr}, with very good agreement.

\begin{figure}
\center{
\includegraphics[width=0.5\columnwidth]{gr.eps}}
\caption{
Scaling of the contact peak of the pair correlation $g(r)$. The theory predicts that in the limit $p\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle\io$ one has
$g(r)/g(1) = \FF(\l)$ with $\l = (r-1) p$. Points are numerical data for $g(r)$ in $d=4$ obtained at several pressures
$p=10^{4},10^5, \cdots, 10^{12}$;
the full line is the theoretical prediction.
}
\label{fig:gr}
\end{figure}

\section*{Discussion}
We have described the marginal phase that is present below
the Gardner transition for hard spheres in $d=\io$.
Using this result we have shown that the jamming transition happens inside the marginal phase and
its criticality in low-dimensional systems
is perfectly described by our approach.
This analysis opens the way for analytically determining many other properties of jammed packings,
such as their shear modulus~\cite{BW06,Yo12},
and the properties of avalanches~\cite{LMW10}.
It also offers a pathway for understanding other glass properties. The nature of aging, for instance, undergoes an abrupt change
at a temperature lower than the glass transition~\cite{St77}, which could signal the crossing between the stable and the fractal glass regimes~\cite{Ri13}.
The observation of dynamical heterogeneities in low temperature glasses~\cite{VKBZ02,BDC08} at timescales much shorter than the inter-basin relaxation
(a fact that does not fit in the two-timescale picture)  may also be taken as a  signature of
the fractal phase, wherein the dynamical correlation length diverges.
Finally, note that in the $d=\infty$ limit, the barriers controlling the inter-basin relaxation should scale as $d$,
and the largest separating subbasins as $d^{1/3}$, suggesting the existence of quasi-localized excitations~\cite{XVLN10}
(stringlike, in three dimensions~\cite{WWW}).

\section*{Methods}

Results are based on the combination of analytical and numerical methods. Analytical results come from the exact solution of hard spheres
in the limit $d\rightarrow} \def\la{\left\langle} \def\ra{\right\rangle\io$, which, for convenience, is obtained
using the replica method, but any other method would give the same result.
The fractal phase is described by a function $\D(y)$ for $y\in [1,y_{\rm max}]$, as in the SK model~\cite{Pa83}. The cutoff $y_{\rm max} \sim p$ diverges with pressure.
With these definitions $\D_{\rm EA} = \D(y_{\rm max})$ is the mean square displacement in the smallest subbasins, where
$\D(y)$ can be computed by numerically solving a set of coupled integro-differential equations obtained from the replica approach (Supplementary Note~1). Numerical results are obtained by standard event-driven
molecular dynamics simulations in $d=3$ to 8~\cite{SDST06,CCPZ12}. Compressions are made using
the Lubachevsky-Stillinger algorithm~\cite{SDST06} (Supplementary Note~2)

\emph{Acknowledgments --}
We thank Silvio Franz for an extremely useful discussion. We also thank Ludovic Berthier, Carolina Brito, Florent Krzakala, Edan Lerner, Federico Ricci-Tersenghi, Tommaso Rizzo, Hajime Yoshino and Matthieu Wyart for many interesting exchanges related to this work. Financial support was provided by
the European Research Council through ERC grant agreement
no. 247328 and ERC grant NPRGGLASS. PC acknowledges support from the Alfred P. Sloan Foundation.


\section{Introduction}
A theoretical wetting study is presented with calculations of the line tension and the boundary tension of a model colloid-polymer mixture. The wetting properties of demixed colloid-polymer mixtures adsorbed at a glass wall have been investigated experimentally and theoretically, following the development of the confocal scanning laser microscopy technique for this system \cite{CSLM1,CSLM2,CSLM3}. This paper is an extension of a previous report in which only the line tension at first-order wetting was investigated \cite{YvesIndekeu}. Now, the line tension well into the partial wetting regime and the boundary tension along prewetting are considered for the same colloid-polymer model.

The line tension is the excess free energy per unit length attributed to the contact line where three immiscible phases meet \cite{Indekeu,Indekeu2}. Here, the three phases are colloidal gas (a phase rich in polymers and poor in colloids), colloidal liquid (rich in colloids and poor in polymers), and a solid (substrate, which is a spectator phase). The solid phase is in mechanical equilibrium with the adsorbates and Young's law is applicable at the three-phase contact line \cite{Bonn}. Experimentally, a wetting transition of what seems to be first-order has been observed for a colloid-polymer mixture with a size ratio $q\approx 1$ \cite{Wijting,Wijting1}. However, to make a definitive statement about the order of this wetting transition more evidence is needed. Also, for size ratio $q\approx 0.6$ only complete wetting states have been observed \cite{Aarts,Aarts1}.

There is a possible extension of the first-order wetting transition off-coexistence where distinct surface phases are present \cite{Indekeu,Indekeu2,Perkovic,Boundarytension}. Here, only the colloidal gas is a stable bulk phase. However, microscopic and mesoscopic films of a colloidal liquid-like phase can be present at the wall. Coexistence of the two thin-film states is possible at a prewetting transition. A line of prewetting transitions starts at the wetting transition and ends in a prewetting critical point. Along this prewetting line, a boundary tension can be defined, which is conceptually related to the line tension. Two possible surface phases share a one-dimensional boundary, and the excess free energy per unit length attributed to this boundary is the boundary tension. In the prewetting critical point the boundary tension vanishes, whereas it takes the same value as the line tension in the wetting transition.

This paper is set up as follows. First, the Cahn-Landau mean-field functional for colloid-polymer mixtures ~\cite{YvesIndekeu,YvesPhD,Inttensionwetting} is reviewed briefly and the interaction potentials and parameters are explained. Then, the interface potentials for the system at first-order wetting, in the partial wetting regime and along prewetting are presented and attention is given to the choice of the Fisher-Jin crossing constraint. The interface potentials are used to calculate the line and boundary tensions. Finally the results are summarised.

\section{A colloid-polymer mixture near a hard wall}
The colloid-polymer mixture near a hard wall can be described by the surface free energy functional
\begin{eqnarray}\label{CPMnearwall}
\gamma[\rho] = \int_{0}^{\infty} \mathrm{d} z \, \Bigl( f(\rho) -
\mu_c \rho(z) + p_c + m(\rho) \bigl( \frac{\mbox{d} \rho}{\mbox{d}
z} \bigr)^2 \Bigr) \nonumber \\ - h_1\rho_1 - \frac{1}{2}g
{\rho_1}^2.  \hspace{4cm}
\end{eqnarray}
The order parameter, $\rho(z)$, is the mean-field colloid number density at perpendicular distance $z$ from the wall, which is located at $z=-\sigma_c/2$, with $\sigma_c$ the colloid particle diameter. This means that the closest position the colloid can have is $z=0$. The colloid density is related to the dimensionless volume fraction through $\phi_c=(\pi/6)\sigma_c^3\rho$. The free energy functional $\gamma[\rho]$ is expressed in the Cahn-Nakanishi-Fisher form and consist of three parts. The first part is composed of the first three terms of the integrand, $f(\rho)-\mu_c\rho(z)+p_c$, which represent the negative of the excess pressure. In this relation, $\mu_c$ and $p_c$ are the equilibrium chemical potential and pressure at two-phase coexistence, respectively. The colloid-polymer mixture is treated as colloidal hard spheres and penetrable polymer spheres. Therefore, the free volume theory with ideal polymers without curvature effects can be applied. The free energy can be written as ~\cite{Phasstruct1}
\begin{equation}\label{Freevolfundamental}
F(N_c,V,T,z_p) = F_0(N_c,V,T) - p^R \alpha(\phi_c) V.
\end{equation}
Here $F_0(N_c,V,T)$ is the hard-sphere free energy, $p^R$ is the polymer reservoir pressure of ideal polymers -- which can be expressed in terms of the reservoir polymer volume fraction $\phi^r_p$ by $p^R=(6/\pi)\sigma_p^{-3}\phi_p^rk_BT$, with $\sigma_p$ the polymer ``diameter" related to the radius of gyration $R_g=\sigma_p/2$ -- and $\alpha(\phi_c)$ is the free-volume fraction available to polymers, which is a function of the colloid volume fraction, $\phi_c$. For $\alpha(\phi_c)$ an approximate expression is found in scaled particle theory ~\cite{Lekkerkerker,Lebowitz}. The hard-sphere free energy is given by the well known Carnahan-Starling expression \cite{CarnahanStarling}.

The second part of the free energy functional is the squared gradient term in the integrand, $m(\rho)(\frac{d\rho}{dz})^2$. The squared gradient is the leading term in the expansion of inhomogeneities of the order parameter \cite{Evansliqvap}. The coefficient of the squared gradient is given by $  m(\rho)\beta = \frac{\pi}{3} \int_{0}^{\infty} dr \,r^4 c(r,\rho) $ where $\beta=\dfrac{1}{k_BT}$ and $ c(r,\rho) $ represents the direct Ornstein-Zernike correlation function with colloid center-to-center distance $r$ and density $\rho$. The direct
correlation function is approximated by \cite{Theoryliquids,BradEvans}
\begin{equation}
c(r, \rho ) = \left\{
    \begin{array}{ll}
        0, & r \leq \sigma \\
        -\beta u(r), & r > \sigma,
    \end{array}
\right. \hspace{1.9cm}
\end{equation}
where $u(r)$ is an attractive interaction potential and arises from the overlap of the depletion zones around the colloid where the polymer is excluded. The two-colloid-particle interaction potential is given by the Asakura-Oosawa depletion potential \cite{AOAO,AOVrij}.

The terms outside the integrand, $-h_1\rho_1-\frac{1}{2}g\rho_1^2$, form the ``wall'' part of the surface free energy functional. These terms correspond to the contact interaction of the colloid-polymer mixture with the hard wall. Here, $h_1$ is the surface excess chemical potential, $g$ is the surface enhancement parameter, and $\rho_1=\rho(0)$, the colloid density at the wall. In the framework of the free-volume theory, $h_1=-\int_{0}^{\infty}dz\;U_2(z)$, where $U_2$ $=p^RV_2$, with  $V_2$ the volume overlap of the depletion zones of a colloidal particle and the hard wall \cite{Inttensionwetting,Aarts}. Also, $g=\int_{0}^{\infty}dz\int{d\vec{r}\;U_3(z,\vec{r})}$, where $U_3$ $=p^RV_3$, with $V_3$ the triple overlap of the excluded volumes of two colloids and the hard wall \cite{strictly}. Recently, Blokhuis and Kuipers argued that this representation of the surface enhancement parameter is incomplete \cite{Blokhuis}. They showed that an extra term has to be included, related to the missing colloid-colloid interactions near the wall. This has a significant impact on the physics of the colloid-polymer mixture. Blokhuis and Kuipers showed that inclusion of this term in the surface enhancement parameter causes wetting reversal at certain values for the aspect ratio $q=\dfrac{\sigma_p}{\sigma_c}$. This means that a drying transition would be observed instead of a wetting transition. Also, the nature of the wetting transition changes as it becomes a second-order phase transition, which leads to the absence of a prewetting regime and consequently no boundary tension can be defined. Another consequence of a second-order wetting transition would be that the line tension at the transition is zero. Although the inclusion of the extra term may constitute a more complete representation of the colloid-polymer mixture, the predictions exhibit poor agreement with experiments \cite{Blokhuis}. For this reason, and in view of our goal to complete the previous account \cite{YvesIndekeu} in which no consideration was made of this extra term in the surface enhancement parameter, this refinement to the theory is ignored. This means that a first-order wetting transition is considered here with a prewetting regime. Thus, for this study the free volume theory with ideal polymers as presented by Aarts et al. \cite{Inttensionwetting} is adopted. A phase diagram of this colloid-polymer system, including a prewetting line, was already presented there \cite{Inttensionwetting}. A thorough generic phase diagram for first-order wetting with representations of all the possible surface phases has been presented by Perkovi\`{c} et al. \cite{Boundarytension}.

\section{Interface Potentials}

The Fisher-Jin interface potential \cite{FisherJin} is the excess free energy per unit area for a homogeneous liquid film of fixed thickness $\ell$ adsorbed at the spectator phase. Therefore, it can be regarded as a constrained, non-equilibrium surface free energy. The interface potential can be obtained from computing the surface free energy $\gamma[\rho]$ for the constrained density profile, the constraint being a crossing density $\rho^{\times}$ at a distance $\ell$ from the surface. An example of a density profile obeying the Fisher-Jin crossing constraint is shown in Fig. \ref{fig:crossing}.
\begin{figure}[h]
  \centering
  \includegraphics[width=0.45\textwidth]{Fig1.eps}
  \caption{An example of a density profile on which the crossing constraint has been imposed. The profile starts with density $\rho_1$ at the wall and is constrained up to the crossing density $\rho^{\times}$. From there the constant of motion is set to zero and the density profile is the unconstrained one, ultimately reaching the bulk vapour phase density $\rho_v$. Note the discontinuity of the first derivative of $\rho(z)$ at $z=\ell$.}
  \label{fig:crossing}
\end{figure}
In the previous report \cite{YvesIndekeu}, it was argued that single-crossing profiles are insufficient to map the complete interface potential. For single-crossing profiles, the density follows the constrained profile up to the value of the crossing density with the corresponding constant of motion. After the crossing density is reached the density profile is the unconstrained profile, for which the constant of motion is zero. The double-crossing profiles follow the unconstrained equilibrium phase portrait only after the crossing density is encountered for the second time. For a graphical explanation of the single-crossing and double-crossing profiles, see Fig. 1b in Ref. \cite{YvesIndekeu}.

Here, we argue that the double-crossing profiles can be avoided if the crossing density is chosen more judiciously. In this section the interface potentials for partial wetting and for prewetting are considered.
\subsection{Partial Wetting}

The interface potential of a colloid-polymer mixture with $q=1$ near first-order wetting with $\phi_c^{\times}=0.03$ -- where the crossing volume fraction is related to the crossing density through $\phi_c^{\times}=(\pi/6)\sigma_c^3\rho^{\times}$ -- was already reported in Ref \cite{YvesIndekeu,YvesPhD}. The necessary equations are repeated briefly. For the optimal density profile the Euler-Lagrange equation has to be solved, while complying with the boundary conditions and crossing constraints. This means that after the density profile reaches the crossing density the constant of motion is zero, whereas leading up to the crossing density the constant of motion can be either positive or negative.
The Euler-Lagrange equation reads
\begin{equation}
2m\dfrac{d^2 \rho}{dz^2}=\dfrac{df(\rho)}{d\rho}-\mu_c
\end{equation}
and the following boundary condition applies
\begin{equation}
-h_1-g\rho_1=2m\dfrac{d\rho}{dz}(0).
\end{equation}
The constant of motion is
\begin{equation}
m(d\rho/dz)^2-f(\rho)+\mu_c\rho-p_c=E.
\end{equation}
\begin{figure}[h]
  \centering
    \includegraphics[width=0.45\textwidth]{Fig2a.eps}
    \includegraphics[width=0.45\textwidth]{Fig2b.eps}
  \caption{Interface potentials, in units of $\frac{k_BT}{(\pi \sigma_c/6)^2}$, near the wetting transition ($S=0$), which can be obtained via a graphical construction of the phase portrait \cite{DeGennes,Cahn}. (a) For aspect ratio $q=1$, the interface potentials are incomplete for crossing colloid volume fraction $\phi_c^{\times}=0.03$, where $\phi_c^{\times}=(\pi/6)\sigma_c^3\rho^{\times}$. The interface potentials are complete for $\phi_c^{\times}=0.05$ and higher. (b) For aspect ratio $q=0.6$, the interface potentials become complete for crossing volume fraction $\phi_c^{\times}=0.09$ and higher. The arrows indicate the evolution of the control parameter $\phi_c^{\times}$ from low to high.}
\label{fig:FOwet}
\end{figure}
\\
Double-crossing profiles were necessary for $\phi_c^{\times}=0.03$ and $q=1$, as the interface potential was incomplete \cite{incomplete}. In Fig. \ref{fig:FOwet} the interface potentials of single-crossing density profiles at the wetting transition, where the so-called spreading coefficient $S\equiv \gamma_{SG}-(\gamma_{SL}+\gamma_{LG})$ is zero, are plotted. It can be seen in Fig. \ref{fig:FOwet}(a) that raising $\phi_c^{\times}$ to $0.05$ takes care of the incomplete interface potential and allows one to obtain a fully smooth and complete one. Further increasing the crossing density causes a shift of the interface potential to lower values of $\ell$. However, the shape of the interface potentials remains the same. The local minimum associated with the thin film starts out at a positive $\ell$ for low values of the crossing density. As the crossing density increases the location of the minimum shifts towards lower values of $\ell$, even reaching negative values. This is not problematic, since the collective coordinate $\ell$ is not the physical thickness of the wetting layer, but a mathematical location which may lie behind the wall \cite{BonnJosPos}. In such cases, in the strictest sense, there is no wetting film, but rather just some modest adsorption at the wall. For first-order wetting and $q=0.6$, which is depicted in Fig. \ref{fig:FOwet}(b), the interface potentials become completely smooth functions for $\phi_c^{\times} \approx 0.09$. High crossing densities lead to familiar interface potentials \cite{BonnJosPos,Indekeu}, with an ordinary barrier without a discontinuity at the maximum between the thin-film minimum and the complete-wetting layer. For low crossing densities, a second branch near the local maximum (at $\ell/\sigma_c \approx 1$) is visible. These branches should be ignored because only the minimal value of $V$ is of interest. The incomplete interface potentials are disregarded for the calculation of the line tension. In this respect, our present approach is a refinement of that taken in \cite{YvesIndekeu}.\\

\begin{figure}[h]
\centering
\includegraphics[width=0.45\textwidth]{Fig3a.eps}
\includegraphics[width=0.45\textwidth]{Fig3b.eps}
\caption{Overview of interface potentials, in units of $\frac{k_BT}{(\pi \sigma_c/6)^2}$, in the partial wetting regime for (a) $q=1$ at $\phi_c^{\times}=0.08$  and (b) $q=0.6$ at $\phi_c^{\times}=0.13$. The arrows indicate the evolution of the control parameter $\phi_p^r$ from low to high.}
\label{fig:InterfaceWetting}
\end{figure}

The interface potentials for states in the partial wetting regime are plotted in Fig. \ref{fig:InterfaceWetting}. For the Fisher-Jin crossing volume fraction values $\phi_c^{\times} = 0.13$ and $0.08$ for $q=0.6$ and $q=1.0$, respectively, the following fitting curve is proposed as an instructive analytical representation of the interface potential.

\begin{eqnarray}\label{VlfitH}
V(\ell) =  \left\{
    \begin{array}{lll}
        a_2(\ell - \ell_1)^2 + a_3(\ell - \ell_1 )^3 + \\
        a_4(\ell-\ell_1)^4, &   0 < \ell < \ell_a \\ \\
        b_1 e^{-\ell/\xi} + b_2 e^{-2\ell/\xi}
         + b_0, & \ell_a < \ell
    \end{array}
\right. \nonumber \\
\end{eqnarray}
This fit consists of two parts. The first part is a Taylor expansion about the thin-film minimum at $\ell_1$ up to third or fourth order. The second part reflects the exponential tail converging to minus the proper spreading coefficient, $-S=b_0$. This tail is approximated by $e^{-\ell/\xi}, e^{-2\ell/\xi}, \cdots$ with $\xi$ the bulk correlation length. This fit introduces an artificial singularity at $\ell=\ell_a$, whereas the interface potentials are smooth throughout. Table \ref{Interfacepotfit1} reports on some typical fitting parameter values near and far from the wetting transition.

\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
\textbf{$\phi_p^r$} & \textbf{$q$} & \textbf{$\ell_{1}/\sigma_c$} & \textbf{$a_2$}   & \textbf{$a_3$} & \textbf{$a_4$} & \textbf{$\ell_{a}/\sigma_c$} \\
\hline
1.917 &  1    & -0.19718  & 0.15357 & -0.08267 & 0.00191 & 1.0044 \\
\hline
3.0  &  1  &  -0.2611   & 0.40599 & -0.19963 & 0.01228 & 1.0603 \\
\hline
\end{tabular}
\end{center}

\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
\textbf{$\phi_p^r$} & \textbf{$q$} &  \textbf{$b_0$}   &   \textbf{$b_1$}   & \textbf{$b_2$}    & \textbf{$\xi/\sigma_c$}    \\
\hline
1.917 &  1    &  0.000041 & 0.77398  & -1.794  & 0.70779  \\
\hline
3.0  &  1  &    0.18356   &  1.40738  & -4.18225   & 0.7226 \\
\hline
\end{tabular}
\end{center}
\caption{Illustrative $V(\ell)$-fitting parameters in the partial wetting regime for $q=1$ ($\phi_c^{\times} = 0.08$).}\label{Interfacepotfit1}
\end{table}

\subsection{Prewetting towards complete wetting}

The interface potential that describes the prewetting transition takes the form
\begin{equation}
V_{pw}(\ell) = V(\ell) + \gamma_{LG}H \ell/\xi,
\end{equation}
where $V(\ell)$ is, with minor modifications, the interface potential at the wetting transition and $H$ is a dimensionless field which takes the system off of two-phase equilibrium and $\gamma_{LG}$ is the liquid-gas surface tension at zero field. The factor $\gamma_{LG}H/\xi$ takes the value of the constant of motion at large $\ell$. Along the prewetting line as $H \rightarrow 0$ the thick-film thickness increases to infinity, and the interface potential converges to the first-order wetting interface potential $V(\ell)$. As can be seen in Fig. \ref{fig:PreWet} the interface potential in the prewetting regime displays two minima, congruent with two surface phases. The first minimum corresponds to the thin film and the second corresponds to the thick film. Table \ref{Interfacepotfit2} reports typical fitting parameter values for several state parameters along prewetting. Note that the Taylor expansion around the minimum is truncated after third order.\\
\\

\begin{figure}[h]
\centering
\includegraphics[width=0.45\textwidth]{Fig4a.eps}
\includegraphics[width=0.45\textwidth]{Fig4b.eps}
\caption{Overview of interface potentials, in units of $\frac{k_BT}{(\pi \sigma_c/6)^2}$, along prewetting for (a) $q=1$ at $\phi_c^{\times}=0.08$  and (b) $q=0.6$ at $\phi_c^{\times}=0.13$. The arrows indicate the evolution of the control parameter $\phi_p^r$ from low to high.}
\label{fig:PreWet}
\end{figure}

\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\textbf{$\phi_p^r$} & \textbf{$q$} & \textbf{$\ell_{1}/\sigma_c$} & \textbf{$a_2$}   & \textbf{$a_3$} & \textbf{$\ell_{a}/\sigma_c$} \\
\hline
1.2 &   1   & 0.081   & 0.0319 & -0.02567 &  0.9087 \\
\hline
1.91 &  1    & -0.19623  & 0.15034 & -0.0782 &   1.1538\\
\hline
\end{tabular}
\end{center}

\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline
\textbf{$\phi_p^r$} & \textbf{$q$} &  \textbf{$b_0$}   &   \textbf{$b_1$}   & \textbf{$b_2$}    & \textbf{$\xi/\sigma_c$} & \textbf{$\gamma_{LG}H/\xi$}   \\
\hline
1.2 &  1 &  -0.08969   &  0.35502  & -0.44494 & 0.81564 & 0.03125 \\
\hline
1.91 & 1  &   -0.001  &  0.78637   & -1.86227  & 0.70257 & 0.0001429 \\
\hline
\end{tabular}
\end{center}
\caption{Illustrative $V_{pw}(\ell)$-fitting parameters along the prewetting line for $q=1$ ($\phi_c^{\times} = 0.08$). $\gamma_{LG}H/\xi$ is given in terms of $\frac{k_BT}{(\pi/6)^2\sigma_c^3}$.}\label{Interfacepotfit2}
\end{table}

\section{Line and boundary tension}

It was previously explained that the excess free energy per unit length of the surface inhomogeneity can be treated within the interface displacement theory \cite{YvesIndekeu}. The line tension is a functional of the displacement of the interface perpendicular to the substrate, and minimising this functional yields the equilibrium line tension \cite{Indekeu1,Indekeu2}. In the squared gradient approximation, this leads to the following expression for the line tension \cite{Indekeu2}
\begin{equation}
\tau=(2\gamma_{LG})^{1/2}\xi \int_{\ell_1}^{\infty}\mathrm{d}\tilde{\ell}\Big[V(\tilde{\ell})^{1/2}-(-S)^{1/2}\Big],
\label{eq:linetension}
\end{equation}
where $\gamma_{LG}$ is the liquid-gas surface tension, $\tilde{\ell}=\ell/\xi$ with $\xi$ the correlation length, and $-S$ is minus the spreading coefficient, which takes the value of $\lim_{\ell\rightarrow \infty}V(\ell)$. The lower integration limit, $\ell_1$, is the thin-film minimum.
\begin{figure}[h]
\centering
\includegraphics[width=0.45\textwidth]{Fig5.eps}
\caption{Line tension, in units of $\frac{k_BT}{(\pi /6)^2\sigma_c}$, as a function of the crossing volume fraction for contact angles $\theta$ (degrees) up to about $40\degree$ for both $q=1$ (full line) and $q=0.6$ (dotted line).}
\label{fig:Linetension}
\end{figure}

Via the spreading coefficient the contact angle can be calculated. In Fig. \ref{fig:Linetension} the line tension is plotted versus the Fisher-Jin crossing volume fraction for several values of the contact angle. The line tension is only weakly dependent on the chosen constraint, as was expected from the interface potentials seen in Fig. \ref{fig:FOwet}. It can be observed that for $q=1$ the line tension rises slightly with increasing values of the crossing density. Farther in the partial wetting regime, as the contact angle increases, the dependence of the line tension on the crossing density becomes more apparent. For $q=0.6$ variation of the line tension with the chosen crossing density is almost non-existent, as it should be in a reliable theoretical model, as it almost looks like a flat line regardless of the proximity to the wetting transition.

For typical colloidal diameters of $10-100 nm$, the computed line tensions are of the order $10^{-14}-10^{-12} N$ near first-order wetting, pass through zero and reach minus their values at first-order wetting for larger contact angles. The line tension values for $q=0.6$ are lower than those for $q=1$, especially near the wetting transition. This can be easily comprehended from the larger area under the $V(\ell)$ curve, as already pointed out near first-order wetting \cite{YvesPhD}. However, with increasing contact angle the difference in line tension decreases and for $\theta\approx 40\degree$ the line tensions are equal. It can be observed in Fig. \ref{fig:Linetension} that the crossing density necessary to provide a complete interface potential, and thus a value for the line tension, increases with increasing contact angle.

\begin{figure}[h]
\centering
\includegraphics[width=0.45\textwidth]{Fig6.eps}
\caption{Boundary ($\tau_B$) and line tension ($\tau$), in units of $\frac{k_BT}{(\pi /6)^2\sigma_c}$, for $q=1$ at $\phi_c^{\times}=0.08$ and $q=0.6$ at $\phi_c^{\times}=0.13$. The boundary tension in the prewetting critical point is zero. The line and boundary tensions are maximal and equal at the wetting transition.}
\label{fig:Boundarytension}
\end{figure}
Like for the line tension, an expression for the boundary tension can be found by using the interface displacement model. This leads to \cite{Indekeu2}
\begin{equation}
\tau_B=(2\gamma_{LG})^{1/2}\xi\int_{\ell_1}^{\ell_2}\mathrm{d}\tilde{\ell}(V_{pw}(\tilde{\ell}))^{1/2}.
\label{eq:boundarytension}
\end{equation}
Here, $\ell_1$ and $\ell_2$ are the thin-film and thick-film minima, respectively. The boundary tension for $q=1$ and $\phi_c^{\times}=0.08$ as well as for $q=0.6$ and $\phi_c^{\times}=0.13$ is plotted against the reservoir volume fraction and the familiar lambda shape \cite{Indekeu,Perkovic,Boundarytension} is found. In the prewetting critical point, the thin and thick-film minimum come together and the boundary tension is zero. The boundary tension reaches the line tension at the wetting transition upon increasing the reservoir volume fraction. {red}The physical parameters at the critical point, the wetting transition, and the prewetting critical point are summarised in Table \ref{tab:summary}.
\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
$q$ & $\phi_p^r$ (cp) & $\phi_c$ (cp) & $\phi_p^r$ (pwcp) & $\phi_c$ (pwcp) & $\phi_p^r$ (wet) & $\tau$\\
\hline
0.6 & 0.4879 & 0.1878 & 0.6007 & 0.0296  & 0.9086 & 0.2134 \\
\hline
1.0 & 0.6364 & 0.1040 & 0.9552 & 0.00681 & 1.917 & 0.7360 \\
\hline
\end{tabular}
\end{center}
\caption{For $q=1$ and $q=0.6$, the location of the critical point (cp), prewetting critical point (pwcp), and wetting transition (wet) are listed in terms of the polymer and colloid volume fractions $\phi_p^r$ and $\phi_c$. Also the value of the line tension $\tau$ at wetting is given, in units of $\frac{k_BT}{(\pi/6)^2\sigma_c}$. }
\label{tab:summary}
\end{table}

\section{Conclusion}

For an adsorbed colloid-polymer mixture, a previous study of the interface potentials and the line tension at the wetting transition using a Cahn-Nakanishi-Fisher functional and an interface displacement model was extended to include the partial wetting and the prewetting regime. The colloid-polymer mixtures are modelled with the free-volume theory applied to ideal polymers without curvature effects. Previously, it was argued that double-crossing profiles were necessary to remedy incomplete interface potentials. Here, it was shown that, near the wetting transition, by increasing the crossing density up to $\phi_c^{\times}=0.05$ and $\phi_c^{\times}=0.09$ for $q=1$ and $q=0.6$, respectively, a smooth and complete interface potential can be obtained. Increasing the crossing density even further brings about a shift of the interface potential to lower $\ell$. Farther in the partial wetting regime, the minimal crossing density to obtain complete interface potentials increases.

At the wetting transition and $q=1$, the line tension increases only slightly with increasing crossing densities. This behaviour becomes more pronounced further into the partial wetting regime, where the line tension becomes negative. For $q=0.6$ the line tension curve is almost constant at the wetting transition and throughout the partial wetting regime. The computed line tensions have a value of $10^{-14}-10^{-12} N$ near first-order wetting, pass through zero and become negative for larger contact angles. For a fixed particle size, close to the wetting transition, the line tension for $q=0.6$ is significantly lower than for $q=1$, whereas they are equal at a contact angle of approximately $40 \degree$. The computed boundary tension along prewetting ranges from zero at the prewetting critical point to the line tension value at the wetting transition.
\end{spacing}

\section{Acknowledgements}
J.K and J.I are supported by KU Leuven Research Grant OT/11/063.

\newpage


\section{Introduction}

Recent experimental research of diffusion and trapping of single particles of $\mu$-
and nm-scales in confined liquid films \cite{Bevan2000,Schuster,Kihm,Banerjee2005,Huang,Sharma} show that
this process presents non-trivial behavior considerably different from a usual
Brownian diffusion in a bulk due to strong interaction of ultrathin fluid matrix
with a base solid interface. The detailed single-molecule tracking
detects the dependence of diffusivity on a distance of layer from the interface.
It should be noted that a majority of evaluated experiments deals with the study of hydrodynamic interaction of colloid
particles with a rigid wall that leads to change of both vertical and horizontal diffusivities \cite{Happel}.
However, the very important parameter of such a system is a particle's size \cite{Bevan2000}. It has been found
\cite{Carbajal,Banerjee2005,Huang,Sharma} that the experiments with micrometer-scale spherical particles diffusing
in a uniform fluid close to the wall quite satisfactory admit the two-dimensional depth-dependent
anisotropic diffusion approach of \cite{Happel}.

At the same time, it was mentioned \cite{Kihm} that particle size change from $\mu$-m to nm scales leads to large deviations from hydrodynamic theory.
Particularly, in the experimental work \cite{Kihm}  it was demonstrated that
observed vertical diffusion coefficient is smaller. The discussion in \cite{Kihm} attributes this fact to the prevalence of intermolecular
forces over hydrodynamical ones for such scales. From this point of view, very important experimental study was evaluated in the work \cite{Schuster},
where extremely small (molecular-size) markers were traced. It has bean detected there that the lateral diffusion is anomalous (superdiffusion) for short times and crosses over
to normal diffusion for long times. This phenomenon was discussed in a connection with the fact that fluids close to wetted interfaces are
know to form a layered structure.

The motion of molecules in such layered structures can be described within the following random walk
picture. The layers close to the interface are practically frozen and
the lateral motion in these layers is slow.
However, particles can diffuse (jump) vertically to layers,
where the  random walks are faster. Thus, the
transient diffusive behavior is based on a coexistence of these
vertical and lateral jumps. Note that the process of molecular diffusion
in a wetting layer differs from the one of macroscopic particles near the wall
which can be described using the hydrodynamical picture \cite{Happel}, {\it vide infra}.

In what follows we adopt a coarse-grained description and consider the
continuous model of layered diffusion: The particle performs diffusive
motion in $y$-direction. The properties of this motion are independent on
the particle's lateral position due to translational symmetry of the fluid
layer in this direction, while the properties of diffusion in $x$-direction
depend on the particle's distance from the surface (i.e. on its $y$-coordinate).
Our model thus corresponds to a one of the layered system.

The models of layered diffusion correspond
to cases when the properties of motion in one ($y$) direction (normal to layers)
do not depend on the coordinate $x$, while the motion within layers (in $x$-direction) depends
on $y$. The properties of motions in $x$ and $y$-directions
differ, defining different models of layered diffusion:
Thus the combination of continuous drift in $y$ with random velocity model in $x$
defines Taylor's model of ``Diffusion by continuous motions'' \cite{Taylor}, the combination
of diffusion in $y$ with the same random velocity model in $x$ defines
a superdiffusive Matheron---De~Marsily model \cite{Matheron}, and combination of diffusion
in $y$ with transport with deterministic mean velocity proportional to $y$
defines the model of diffusion in a shear flow introduced by Novikov \cite{Novikov}.
More general $y$-dependences of $x$-velocity define a class of models in \cite{BenNaim}.

The model of a layered system is not limited to hydrodynamical situations.
Thus, a similar model was recently used
for modeling of superdiffusion in a small world network
\cite{Postnikov}. Namely, the set of slow walks along a regular
background lattice and fast jumps along ``short links'' was replaced
by the set of simple diffusion equations with the hierarchy diffusion
coefficients depending on the distance between background nodes.

\section{The model}

A simple model described here adds an interesting example to the theory
of diffusion in layered systems and gives rise to an effective equation
of the Batchelor's type for lateral diffusion. Due to symmetry in the directions
parallel to the surface, it is enough to consider one lateral coordinate, which
is denoted $x$ in what follows. For this reason we consider 2D picture at first.
Then we discuss the generalization on 3D case, which is quite straightforward.
Our model corresponds thus to a particle performing diffusion in $x$ and $y$-directions,
where the properties of the vertical diffusion (in $y$-direction) do
not depend on $x$. On a coarse-grained level (the steps of a random walk are not
resolved) this model is described by a pair of Langevin equations
\begin{eqnarray}
\dot{y} &=& \sqrt{2D_y} \xi_y(t) \label{motiony} \\
\dot{x} &=& \sqrt{2D_x(y)} \xi_x(t) \label{motionx}
\end{eqnarray}
with $D_y$ and $D_x(y)$ being the diffusion coefficients in $y$- and
in $x$-directions, and $\xi_y(t)$, $\xi_x(t)$ being Gaussian white
noises with zero mean and correlation property given by
$\langle\xi_\alpha(t) \xi_\beta (t') \rangle = \delta_{\alpha\beta}
\delta(t-t')$. Note that according to Eqs. (\ref{motiony}) and
(\ref{motionx}) the motion in $y$ is independent on the motion
in $x$-direction, while the motion in $x$ does depend on the
actial $y$-coordinate through $D_x(y)$ and thus depends on the prehistory of $y$-motion.
The random variables $x$ and $y$ are not independent.

The evolution of two-dimensional joint probability
density $p(x,y,t)$ for marker's position is then described by the
Fokker-Planck equation corresponding to the equation of anisotropic
diffusion:
\begin{equation}
\frac{\partial p}{\partial t}=D_x(y)\frac{\partial^2 p}{\partial x^2}+
D_y\frac{\partial^2 p}{\partial y^2}.
\label{init_eq}
\end{equation}
The Dirac delta function
$p(x,y,0)=\delta(x,y+0)$ is used as an initial condition.  We are
interested only in the time dependence of marginal distribution of
$x$, which is experimentally accessible \cite{Schuster}.
We note that although the situation with molecular tracer discussed here is different from the case of
motion of a massive particle in a quiescent fluid, as discussed e.g. in \cite{Sancho82}, since
a molecular motion is never underdamped. However, the final equation is
exactly of the form obtained in this work, due to the fact that it is
the only form corresponding to thermodynamical equilibrium (in the absence
of external forces) in a stationary state.

Although Eq.(\ref{init_eq}) itself does not have to possess a
solution which factorizes into a product of functions which depend
only on one of the two co-ordinates, the solution of our physical problem
given by Eqs.(\ref{motiony}), (\ref{motionx}) has to factorize:
$p(x,y,t)=p(x,t)p(y,t)$, although $x$ and $y$ are not independent.
The physical reason for this is as follows.

First let us remark, that while the two arguments $x$ and $y$ of
$p(x,y,t)$ are the values of the corresponding random variables $x$
and $y$, its third argument, time, is not a random variable but essentially
an additional condition: $p(x,y,t)=p(x,y|t)$ is a conditional joint
probability of $x$ and $y$ provided their measurement took place at
time $t$. Note now that according to Eq.(\ref{motiony}) the motion in
$y$ direction is independent on $x$, and moreover the
initial condition for this motion in $y$ does not depend on
$x$. Therefore the conditional probability of $y$ given $x$ is
essentially a non-conditional one: $p(y|x,t) = p(y|t)$.  According to
Bayes theorem $p(x,y|t) = p(y|x,t)p(x|t) = p_y(y|t)p_x(x|t)$: the $x$
and $y$-displacements are {\it conditionally independent} \cite{Dawid1979}.
Rewriting it in our initial notation we get
\begin{equation}
p(x,y,t)=p_x(x,t)p_y(y,t).
\label{separation}
\end{equation}
Independently on whether this discussion is convincing or not, one can always prove that
the final solutions obtained by the separation ansatz are indeed the solutions of Eq.(\ref{init_eq}),
i.e. that the product of Eq.(\ref{sol1infy}) and Eq.(\ref{sol1infx}) indeed solves Eq.(\ref{init_eq})
in an unbounded system with $D_x(y)=ky$, and that the product of Eq.(\ref{sol1}) and of a symmetric Gaussian in $x$
with the dispersion given by Eq.(\ref{MSD}) solves it for a layer of a finite thickness, as they indeed do.

Integrating both parts of Eq.(\ref{init_eq}) with respect to $x$
we obtain the equation for $p_y(y,t)$:
\begin{equation}
\frac{\partial p_y}{\partial t}=D_y\frac{\partial^2 p_y}{\partial y^2}.
\label{eq_u_x}
\end{equation}

Inserting Eq.(\ref{separation}) into Eq.(\ref{init_eq}) and integrating
over $y$ we obtain
\begin{equation}
\frac{\partial p_x}{\partial t}=\left[\int_0^H D(y) p_y(y,t) dy \right]
\frac{\partial^2 p_x}{\partial x^2},
\label{eq_u_y}
\end{equation}
with $H$ being the film's thickness,
i.e. the one-dimensional Batchelor's equation \cite{Batchelor}
\begin{equation}
\frac{\partial p_x}{\partial t}=D(t)\frac{\partial^2 p_x}{\partial x^2},
\label{Batchelor}
\end{equation}
with time-dependent diffusion coefficient
\begin{equation}
D(t)=\int_0^HD_x(y)p_y(y,t)dy.
\label{D}
\end{equation}
Note that this Batchelor's equation arises here not as the mean-field type
approximation, but as an exact equation for the corresponding marginal probability
density. Introducing a new variable $\tau=\int_0^tD(t)dt$ allows to reduce Eq. (\ref{Batchelor})
to the simple diffusion equation
\begin{equation}
\frac{\partial p_x}{\partial \tau}= \frac{\partial^2 p_x}{\partial x^2}.
\label{px}
\end{equation}
Note that the change of variables from $t$ to $\tau$ is always possible within the model
since the effective diffusion coefficient $D(t)$ given by Eq.(8) is always non-negative
and therefore its integral over time is a monotonically growing (and therefore invertible) function.
From this it follows
that the corresponding probability density is Gaussian with a
mean square displacement (MSD) $<x^2>=2\tau$.

In a 3D system uniform in $(x,z)$ coordinates with $D_z(y)=D_x(y)$ still depending on the vertical
$y$ position, the system of Eqs.~(\ref{motiony})--(\ref{motionx}) is supplemented by
\begin{equation}
\dot{z}=\sqrt{2D_x(y)} \xi_z(t). \label{motionz}
\end{equation}

Thus, from the Eqs. (\ref{motiony})--(\ref{motionx}), (\ref{motionz}) and the symmetry in the directions
parallel to the surface, it follows that the three-dimensional joint probability
density $p(x,y,z,t)$ satisfies the equation:
\begin{equation}
\frac{\partial p}{\partial t}=D_x(y)\left[\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial z^2}\right]+
D_y\frac{\partial^2 p}{\partial y^2}.
\label{init_eq3D}
\end{equation}

As before, all displacements are conditionally independent
and
the three-dimensional joint probability
density will be factorized as
\begin{equation}
p(x,y,z,t)=p_x(x,t)p_y(y,t)p_z(z,t).
\label{3Dfactorization}
\end{equation}

The equation for $p_z(z,t)$ has exactly the form of the Eq.~(\ref{px}).

Finally, since the motions in both lateral directions are independent, the resulting MSD is
\begin{equation}
<r^2>=<x^2>+<z^2>=2<x^2>=4\tau.
\label{r2}
\end{equation}
In other words, it is simply two times larger then in 2D case, and the definition of $\tau$ follows from the Eq.~(\ref{D}) as above.

\section{Solutions}

\subsection{Superdiffusion in a thick film}

At first, we consider an idealized situation: the half-plane corresponding
to the film of infinite thickness $H=\infty$. In this case the solution of
Eq. (\ref{eq_u_x}) is the Gaussian
\begin{equation}
p_y(y,t)=\frac{1}{\sqrt{D_y\pi t}} \exp \left(-\frac{y^2}{4D_yt}\right).
\label{sol1infy}
\end{equation}
Thereafter, one needs to determine the vertical dependence of the
horizontal component of the diffusion coefficient. As it will be
shown below, the experimental data are well-fitted
by a linear dependence $D_x(y)=ky$. Under this condition, the substation
(\ref{sol1infy}) into (\ref{D}) gives
\begin{equation}
D(t)=2k\sqrt{D_yt/\pi}.
\label{Dtinf}
\end{equation}
Therefore, the averaged lateral mean-square displacement corresponds to
the superdiffusion:
\begin{equation}
\left<x^2\right>=\frac{8k}{3}\sqrt{\frac{D_y}{\pi}}t^{\frac{3}{2}},
\label{MSDinfthick}
\end{equation}
and the corresponding distribution of lateral positions is given by a
Gaussian
\begin{equation}
p_x(x,t)=\frac{\exp \left(-\frac{x^2}{\sqrt{(4/9\pi)k^2D_y^3t^3}} \right)}{\sqrt[4]{(64/9)\pi k^2D_y^3t^3}}.
\label{sol1infx}
\end{equation}

 As it was discussed above, in 3D case additional multiplier $p_z(z,t)$ in the Eq.~(\ref{3Dfactorization})
will have exactly the same form as (\ref{sol1infx}) with replacing $x \to z$, and MSD $<r^2>$ is given by
twice the expression (\ref{MSDinfthick}).

\subsection{Film of finite thickness}

Now let us consider the realistic situation, when the thickness $H$ is
finite. In such a case the equation (\ref{eq_u_x}) with impenetrable
boundaries has a standard solution in the form of the eigenfunction expansion:
\begin{equation}
p_y(y,t)=\frac{1}{H}+\frac{2}{H}\sum_{n=1}^{\infty}
e^{-\frac{D_y\pi^2n^2}{H^2}t}\cos\left(\frac{\pi n}{H}y\right).
\label{sol1}
\end{equation}
Under the same condition $D_x(y)=ky$, the time dependence of the
averaged coefficient of lateral diffusion is represented by the
 series
\begin{equation}
D(t)=\frac{kH}{2}\left(
1-\frac{8}{\pi^2}\sum_{m=0}^{\infty}\frac{e^{-\frac{D_y\pi^2(2m+1)^2}{H^2}t}}
{(2m+1)^2}
\right).
\label{Dt}
\end{equation}
From Eq. (\ref{Dt}) it follows that the effective diffusion
coefficient is initially equal to zero due to sticking of a
marker to the solid substrate, monotonically grows with time and then saturates.
Its final value $D(\infty)=kH/2$ is achieved when the markers redistribute homogeneously
over the whole liquid layer.

Integrating Eq. (\ref{Dt}) over time we
get the desired expression for a mean-square-displacement:
\begin{equation}
<\hspace{-1mm}x^2\hspace{-1mm}>=kH\left(
t+\frac{8H^2}{D_y\pi^4}\sum_{m=0}^{\infty}\frac{e^{-\frac{D_y\pi^2(2m+1)^2}{H^2}t}-1}
{(2m+1)^4}
\right).
\label{MSD}
\end{equation}

Analyzing Eq. (\ref{MSD}), one can see that non-linear growth at a small
time changes to the asymptotic linear one  corresponding to a normal diffusion:
\begin{equation}
<x^2\hspace{-1mm}>_{\infty}=-\frac{kH^3}{12D_y}+kHt.
\label{x2norm}
\end{equation}

Fig.~\ref{MSDfig} shows the typical behavior of the mean squared displacement
as a function of time. It clearly shows the non-stationary character of
diffusion in ultrathin liquid film. During the first stage the
motion is superdiffusive: the curve is well fitted
by power-law function with the exponent $3/2$.
The microscopic physical background of this can be explained
by jumps of markers into higher layers and fast transport there.
The further saturation of layers by markers leads to decrease of the
relative contribution of processes ``jump to a higher layer and fast
diffusion there''. Correspondingly, the growth of the
mean squared displacement decelerates. Finally, the uniform
redistribution of markers changes a kind of transport to a normal
diffusion.

All corollaries concerning simple 3D generalizations are the same as in the previous subsection.

\begin{figure}
\includegraphics[width=\columnwidth]{MSD}
\caption{The time dependence of mean-square-displacement (solid line)
for the parameters: $H=17 \; \mathrm{nm}$, $k=0.4384 \; \mathrm{cm^2/(s\cdot
nm)}$, $D_y=8\cdot 10^{-13} \; \mathrm{cm^2/s}$. Two additional referent
lines correspond to the superdiffusion $<x^2>=0.02t^{3/2}
10^{-10} \; \mathrm{cm^2}$ (dashed line) and the asymptotic normal
diffusion $<x^2>_{\infty}$. (dash-dotted line)}
\label{MSDfig}
\end{figure}

\section{Comparison with experiment}

To show that the model presented does make sense, the comparison with direct measurements
of inhomogeneous molecular diffusion in thin liquid film deposited on
solid surface \cite{Schuster} is carried out. Let us first show, that
the assumption of the linear dependence of the microscopic lateral diffusion
coefficient on the distance from solid substrate holds. The least-square linear fit
of  the experimental data for lateral diffusion coefficients taken from \cite{Schuster}
is presented Fig.~\ref{fit}.  This fit  justifies the assumption the approximate law for
$D_x(y)$ with a rather good accuracy.

\begin{figure}
\includegraphics[width=\columnwidth,height=0.6\columnwidth]{verticalD}
\caption{The linear fit of experimental data of Ref. \cite{Schuster} (Fig. 6b there; black squares corresponds to the same symbols in the cited source), which
confirms the chosen dependence $D_x(y)=ky$ and allows to determine the
tangent $k=0.4384 \; \mathrm{cm^2/(s\cdot nm)}$.}
\label{fit}
\end{figure}

Using the value of the slope obtained, the calculation according to the
formula (\ref{MSD}) are completed for the experimental film thickness
(note that the series (\ref{MSD}) has a relatively fast convergence, e.g. the maximal relative
difference between result presented by truncations at 11th and 10th terms is 2\% and at between
truncated 21th and 20th terms is 0.2\%). Note, that the vertical component of diffusion
coefficient is unknown from the direct measurements. Therefore, it is
estimated in such a way that the resulting curve fits the experimental
points, see Fig.~\ref{MSDexp}. These experimental points represented as stars in Fig.~\ref{MSDexp}.
 Note also that if one need to take into account 3D picture with two-dimensional lateral layers, the curves keep the same shape, but the adjustable parameter would be simple $\sqrt{2}$ times smaller. This follows from the Eqs.~(\ref{D}), (\ref{r2}).

The experimental points correspond to a single trajectory from the set of the overall 4 trajectories
kindly provided by authors of Ref. \cite{Schuster}.
For all of them the behavior at short times (governed by practically deterministic process
of the departure from the surface) is practically the same.
For longer observation times, the individual trajectories show various fluctuation features
like large excursions and returns, and the averaging over our very restricted sample does not produce
better results than taking an individual one. We use the trajectory, MSD of which corresponds
to most regular curve over the time interval considered.

The considered time range allows to cover two specific regimes:
i) During the interval from 0 to approximately 400~ms, the process is superdiffusive
MSD ($<x^2>\sim t^{3/2}$). It coincides with the analytically exact superdiffusion
obtained for semi-infinite region. The natural reason for this is that at short times a walker can
not approach the upper boundary of film. Thus there is no sufficient difference
between solutions as well as for physical experimental picture.  ii) The transient
process from superdiffusion to normal diffusion starting from time approximately 400~ms.
Within this subinterval, the MSD curve can not be represented as a power-law function,
compare with the Fig.~\ref{MSDfig}, however the series solution (\ref{MSD}) adequately
reproduces experimental data.

The asymptotic large-time MSD corresponding to the normal diffusion,
see Fig.~\ref{MSDfig}  for simulated curve and Fig.~\ref{MSDnormdiff} confirming this conclusion by the experimental data. See also discussion of long observation times in \cite{Schuster}.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{MSDexpTheor}
\caption{The comparison of calculated curve of a
mean-square-displacement  (solid line) with experimental data (asterisks) provided by authors of
\cite{Schuster}  for short observation time. Parameters of the film: $H=17 \; \mathrm{nm}$, $k=0.4384
\; \mathrm{cm^2/(s\cdot nm)}$; vertical diffusion coefficients
$D_y=8\cdot 10^{-13} \; \mathrm{cm^2/s}$  (or $8\cdot 10^{-13}/\sqrt{2} \; \mathrm{cm^2/s}$) for 3D case.}
\label{MSDexp}
\end{figure}

 Thus, that anomalous diffusive
behavior in lateral direction is caused by very weak vertical
diffusion, i.e. rare jumps between layers: their rate is five orders of magnitude
smaller then the lateral one.

Note also that our model deals with averaged
motion. Correspondingly, it does not take into account fluctuations
along individual trajectories, which are significant for large times of observation
mirrored in growing fluctuations in $D$ \cite{Schuster}.
Since the total number of trajectories obtained in the experiment is limited, and these trajectories show considerable differences between different runs, as they should. Fitting the results of a theory relying on ensemble-averaged quantities to this data is therefore connected with large statistical uncertainties, and one can mostly speak about trends, which are correctly captured by our theoretical description. The more convincing proof can be only given by performing more experiments in longer runs.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{norm_diff}
\caption{ The comparison of calculated curve of a
mean-square-displacement ($<x^2\hspace{-1mm}>_{\infty}\sim kHt$, see the Eq.~(\ref{x2norm}). The constant vertical shift does not coincide with the analytical expression due to certain irregular markers' jumps and stops during the experiment) with experimental data (asterisks) provided by authors of
\cite{Schuster}  for long observation time. Parameters of the film: $H=17 \; \mathrm{nm}$, $k=0.4384
\; \mathrm{cm^2/(s\cdot nm)}$; vertical diffusion coefficient
$D_y=8\cdot 10^{-13} \; \mathrm{cm^2/s}$ (or $8\cdot 10^{-13}/\sqrt{2} \; \mathrm{cm^2/s}$ for 3D case).}
\label{MSDnormdiff}
\end{figure}

\section{Additional discussion}

One of our basic simplifications is based on the assumption of a
constant vertical diffusion coefficient coexisting with
height-dependent lateral one.  However, detailed studies
of a near-wall motion of mesoscopic particles show that both
components of the diffusion coefficient should depend on the
distance from a wall, see the theoretical derivation and
discussion in \cite{Happel}.
In this case the explicit $y$-dependence of the vertical diffusion coefficient
should be taken into account, and
Eq.~(\ref{init_eq}) should be replaced by a more general one
\begin{equation}
\frac{\partial p}{\partial t}=D_x(y)\frac{\partial^2 p}{\partial x^2}+
\frac{\partial}{\partial y}D_y(y)\frac{\partial p}{\partial y}.
\label{init_eq_cd}
\end{equation}
(If starting on the Langevin level, the kinetic (Hanggi-Klimontovich) interpretation
of the Langevin equation for vertical motion should be used when assuming that the
equilibrium distribution of the molecules in $y$-direction is homogeneous \cite{Sokolov2010}).
Although Eq.~(\ref{init_eq_cd}) is more complex than our initial Eq.~(\ref{init_eq}),
the assumption of independence of motion in $y$ on the lateral position and thus the
conditional independence and the factorization of $x$- and $y$-motions still hold.

Now, let us subdivide the vertical diffusion coefficient in (\ref{init_eq_cd}) into
a constant and a variable parts: $D_y(y)=D_y^{(c)}+D_y^{(v)}(y)$.
Then the Eq.~(\ref{eq_u_x}) takes the form:
\begin{equation}
\frac{\partial p_y}{\partial t}=
\left(D_y^{(c)}+D_y^{(v)}(v)\right)\frac{\partial^2 p_y}{\partial y^2}+
\frac{\partial D_y^{(v)}}{\partial y}\frac{\partial p_y}{\partial y}
\label{eq_u_x_mod}
\end{equation}
Supposing that $D_y^{(v)}<<D_y^{(c)}$ we can neglect the second term in
brackets in the first term in the right hand side of this equation (in a
factor in front of the second derivative) and conclude that the coordinate
dependence of $D_y$ mainly leads to the emergence of a drift
term in $y$-direction. The last one actually changes the motion from a
pure Brownian one to a combination of diffusion and regular drift. These results in a
change of the effective MSD-exponent in $y$ from 1.5 closer to 2 during
the first stage but does not influence the late stage of the process: after a uniform
redistribution of walkers over a bulk the averaged lateral MSD will
correspond to the normal diffusion. The numerical solution of
(\ref{eq_u_x_mod}) with $D_y^{(v)}(y)=0.05\div0.2D_y^{(c)}y$ shows that initial superdiffusion in this
case is sufficiently faster than the experimentally observed one,
and moreover that the transient to a normal diffusion is shorter.
In fact the resulting process comprises short-time regular ballistic drift
($\left<x^2\right>\sim t^2$) leading to the uniform filling of the layer with
succeeding normal lateral diffusion.

Thus, the
assumption of uniform small vertical diffusion coefficient can be considered as
adequate for the description of the observed pattern of  lateral diffusion  for extremely small tracers, confirming the supposition made in \cite{Kihm}. The additional
confirmation of this assumption can be found in the work \cite{Schuster2003},
where the interpretation of experiments is provided. Authors argued
that the microscopic molecular picture of media corresponds to the
layered structure with small probability of jumps of molecules
(nanometer scale objects) between adjacent layers.  Our model
discussion confirms these conclusions.

\section{Summary}
We considered lateral diffusion of a particle in a layered ultrathin fluid medium and obtained the
solutions for the lateral mean-square-displacement. The results demonstrate the transition
from superdiffusion at short times to asymptotically normal diffusion, as
is confirmed by the experimental findings.

From the theoretical point of view, our work discusses an interesting case when the Batchelor's equation in 1D
appears to be exact, emerging from the full two-dimensional Fokker-Plank equation
due to conditional independence of motions along both coordinates. Note that these conditionally independent
motions are still connected via the common parameter: the time.

\section*{Acknowledgments}
We thank Prof.~C.~von~Borzcyskowski and Dr.~J.~Schuster for provided
experimental data and valuable discussions.


\section{Introduction}
\label{Sec:Intro}

A \emph{symmetry} of an object is a geometric transformation which leaves the object unchanged. So, for example, an object with 3-fold rotational symmetry has three symmetries: rotation by $120^\circ$, rotation by $240^\circ$, and the trivial symmetry, where we do nothing. The symmetries of an object naturally form a group under composition.
Care must be taken to differentiate between the symmetry group of an object, consisting of geometric transformations that leave the object unchanged, and the abstract group, which only contains information about how the elements of the group interact with each other under composition.

\begin{figure}[htbp]
\centering

\subfloat[An object with symmetry group isomorphic to $D_4$.]
{
\includegraphics[width=0.3\textwidth]{Figures/diagram_for_224.pdf}
\label{Fig:224}
}
\quad
\subfloat[Another object with symmetry group isomorphic to $D_4$.]
{
\includegraphics[width=0.3\textwidth]{Figures/diagram_for_star44.pdf}
\label{Fig:*44}
}
\quad
\subfloat[A Cayley graph for $D_4$. The edges with arrows correspond to rotations, the other edges correspond to reflections.]
{
\includegraphics[width=0.3\textwidth]{Figures/D4_cayley_graph}
\label{Fig:D4_cayley}
}

\caption{Symmetric designs on the sphere.}
\label{Fig:spherical_symm_exs}
\end{figure}

As an example, consider the objects pictured in Figure \ref{Fig:spherical_symm_exs}. The object shown in Figure \ref{Fig:224} has no planes of mirror symmetry, but has ten \emph{gyration points} (points of rotational symmetry, marked with a small circle).
The symmetry group of this object consists of rotations about these gyration points by multiples of either $90^\circ$ or $180^\circ$ (depending on the kind of gyration point).
The object shown in Figure \ref{Fig:*44} has four planes of mirror symmetry, marked by their intersection with the sphere. The symmetry group of this object consists of reflections in these planes, together with rotations about the gyration points at the top and bottom of the sphere.

These objects have different sets of symmetries, although as abstract groups they are both isomorphic to the dihedral group with eight elements, $D_4 = \langle a,b \mid a^{4} = b^2 = 1, b a b = a^{-1}\rangle$.
Figure \ref{Fig:D4_cayley} shows a Cayley graph which abstractly represents the group $D_4$, although the figure itself does not have $D_4$ symmetry.\footnote{Cayley graphs are purely combinatorial objects with no intrinsic embedding in space associated with them, although they are often drawn symmetrically embedded in space (in this case with order 4 rotational symmetry). Under our definition then, the Cayley graph of a group does not count as a physical object with the symmetry of that group, unless it is embedded in a way that actually has that symmetry.}

As artists, we are particularly interested in
the symmetries of real world physical objects.
Three natural questions arise:
\begin{enumerate}
\item Which groups can be represented as the group of symmetries of some real-world physical object?
\item Which groups \emph{have actually} been represented as the group of symmetries of some real-world physical object?
\item Are there any glaring gaps -- small, beautiful groups that should have a physical representation in a symmetric object but up until now have not?
\end{enumerate}

In one sense, the first question is answered in detail in \emph{The Symmetries of Things} \cite{symmetries_of_things} by Conway, Burgiel and Goodman--Strauss. For example, the spherical symmetries of objects in $\mathbb{R}^3$ with symmetries consisting of rotations and reflections whose axes and reflection planes all pass through a single point are entirely classified.
In another sense things are not so clear -- we can draw pictures of the wallpaper symmetry groups, but real, physical wallpaper never extends infinitely far in the plane. Therefore the design is not actually unchanged by a translation, or even approximately unchanged. However, the viewer can infer how the pattern should continue, and so in this sense the wallpaper groups can be and have been represented by real world objects. Similarly, no sculpture (that we know of) is exactly symmetric at a subatomic level, but the viewer infers the intended symmetry. Using models such as the Poincar\'{e} disk model of hyperbolic space we can make physical representations of hyperbolic symmetry groups, for example in Escher's ``Circle Limit'' series. Likewise, artists have used shadows and projections to represent higher dimensional symmetry groups as physical objects.

To our knowledge, the second question has not been investigated in this generality.\footnote{A related question that \emph{has} been investigated in detail is the issue of which of the 17 wallpaper groups appear in the Alhambra. This is a topic of some debate (see for example \cite{grunbaum_alhambra_symm_groups}), with questions over whether or not different colours can be counted as the same in order to show a symmetry group, how large a design needs to be to represent a group and so on. However, whether or not these symmetry groups appear in the Alhambra, they certainly have appeared in modern works.}
There are of course infinitely many groups, but only finitely many physical objects that have ever been made. So for example, there is some smallest finite N for which there has never been a physical object with N-fold rotational symmetry.
Certainly all of the polyhedral symmetry groups have been represented in sculpture since antiquity, if not represented in nature beforehand. As we move to representations of more esoteric groups (for example requiring a model of hyperbolic space), the number of physical representations surely goes down.

As for the third question, there does seem to be such a gap for the \emph{quaternion group} $Q_8$.\footnote{Note that there are many drawings of the Cayley graph of $Q_8$ with various embeddings, as well as art representing its multiplication table. See for example Gwen Fisher's Quaternion Quilt \cite{fisher}.}
In Section \ref{Sec:Quat group} we give an abstract definition of the quaternion group.
In Section \ref{Sec:Visualise quat group} we give an intuitive way to visualise the action of $Q_8$ on a hypercube, by building the boundary of the hypercube from eight decorated cubes.
In Section \ref{quat group in sculpture} we
give our strategy for producing a sculpture whose symmetry group is $Q_8$, using stereographic projection from the unit sphere in 4-dimensional space. To our knowledge, this is the first sculpture that specifically exhibits $Q_8$ as a symmetry group.
In Section \ref{future directions} we give some directions for future work.

\section{The Quaternion group}
\label{Sec:Quat group}

The quaternion group $Q_8$ is a beautiful group of order eight, not to be confused with the real quaternions or the unit quaternions.
The real quaternions $\mathbb{H}$ are a number system that extends the real numbers, similarly to the way in which the complex numbers also extend the real numbers.\footnote{The symbol $\mathbb{H}$ is chosen in honour of the discoverer of the real quaternions, William Rowan Hamilton.}
$$\mathbb{H} = \{w +xi+yj+zk \mid w,x,y,z\in\mathbb{R}, i^2=j^2=k^2=ijk=-1 \}$$
The quaternion group consists of eight elements of $\mathbb{H}$ under multiplication.
$$Q_8=\{\pm1, \pm i, \pm j, \pm k\}$$

\begin{wraptable}[17]{r}{0.46\textwidth}
\vspace{-15pt}
\centering
\[
\begin{array}{ r | rrrrrrrr }

  & 1 & -1 & i & -i & j & -j & k & -k \\
  \hline
1 & 1 & -1 & i & -i & j & -j & k & -k \\
-1 & -1 & 1 & -i & i & -j & j & -k & k \\
i & i & -i & -1 & 1 & k & -k & -j & j\\
-i & -i & i & 1 & -1 & -k & k & j & -j\\
j & j & -j & -k & k & -1 & 1 & i & -i\\
-j & -j & j & k & -k & 1 & -1 & -i & i\\
k & k & -k & j & -j & -i & i & -1 & 1\\
-k & -k & k & -j & j & i & -i & 1 & -1

\end{array}
\]

\caption{The multiplication table (i.e., Cayley group table) for the quaternion group. The convention for the order of multiplication used here is that the row label is first, followed by the column label.}
\label{mult_table}
\end{wraptable}

We give the multiplication table for $Q_8$ in Table \ref{mult_table}. The quaternion group has one element of order one, one element of order two, and  six elements of order four.
We can view $\mathbb{H}$ as $\mathbb{R}^4$ with an added multiplication structure. Multiplication by a quaternion then acts as a transformation on $\mathbb{R}^4$. The unit quaternions are quaternions $q=w+xi+yj+zk$ with magnitude $|q|=\sqrt{w^2+x^2+y^2+z^2}=1.$ These act by 4-dimensional rotation on $\mathbb{R}^4$, and so $Q_8$ acts by rotation on $\mathbb{R}^4$ (and also acts by rotation on $S^3$, the sphere of radius 1 in $\mathbb{R}^4$). In fact, $Q_8$ also acts nicely on the eight cubical cells of the hypercube, permuting them. To see this, note that the eight cells are given by $\{(1,x,y,z)\mid x,y,z\in\mathbb{R}\}$, $\{(-1,x,y,z)\mid x,y,z\in\mathbb{R}\},$ and similarly but with the $1$ or $-1$ in the three other possible coordinate positions. If, for example, we act on the point $(1,x,y,z) = 1+ix+jy+kz$ on the right by $i$, we get $(1+ix+jy+kz)i = i-x-ky+jz=-x+i+zj-yk=(-x,1,z,-y)$, which is on one of the other cells.

\section{Visualizing the quaternion group}
\label{Sec:Visualise quat group}

\begin{wrapfigure}[14]{l}{0.35\textwidth}
\vspace{-10pt}
\centering

\includegraphics[width=0.35\textwidth]{Figures/monkeyblock_net_and_cube_left_screw.jpg}
\caption{A monkey block.}
\label{monkeyblock_net}
\end{wrapfigure}

In order to visualize the symmetries of the quaternion group, we begin by assembling a hypercube out of eight \emph{monkey blocks}.\footnote{The reason for our use of a monkey will become clear in Section \ref{quat group in sculpture}.} See Figure \ref{monkeyblock_net}. Each monkey block is a cube which, due to the patterns on the faces, has no symmetry itself. One face of the cube has a monkey's face on it with her tongue sticking out to the right, while the opposite face on the cube has a monkey face with her tongue sticking out to the left. In addition, opposite faces are rotated a quarter turn from each other.
A right monkey paw is opposite a quarter-turned left monkey paw, and a question-mark tail is opposite a quarter-turned un-question-mark tail.

Given two of these monkey blocks, one can match any face on one block to its mirrored version on the other. Given an infinite number, one can match up all the right paws with all the left paws to get an infinite line of monkey blocks. See Figure \ref{line}. This infinite line has translational symmetry: the whole line can move forward or backwards any multiple of four blocks. It also has a more unusual symmetry: it can move forward one and rotate $90^\circ$ to the left, giving a screw motion.

\begin{figure}[h]
\centering

\includegraphics[width=0.8\textwidth]{Figures/monkeyblock_line_left_screw_v2.jpg}

\caption{An infinite line of monkey blocks.}
\label{line}
\end{figure}

The next obvious thing to try when matching up cubes might be tiling monkey blocks infinitely in two or three dimensions, but there is no way to match four blocks around a single edge.
However, if we allow ourselves to attach blocks with a $90^\circ$ bend into the fourth dimension, we can take four blocks from our infinite line and attach them in a loop. See Figure \ref{ring}.
In this four-dimensional arrangement (which is not really shown in Figure \ref{ring}!) we have two perpendicular planes of rotation. Rotation in the first plane corresponds to moving the blocks around the circle passing through the paw sides of the monkey blocks in the direction of the arrow (labelled with a ``1''). Rotation in the second plane (which is harder to see in the figur
) corresponds to keeping each block in place but spinning it on the circle (the arrow labelled ``2'').

The symmetry of the ring of blocks consists of simultaneous\footnote{Two perpendicular rotations in four dimensions may be performed in either order to obtain the same result.} $90^\circ$ rotation in the two perpendicular planes, something not possible in fewer than four dimensions. This makes the symmetry an \emph{isoclinic rotation}: two equal and perpendicular rotations. Isoclinic rotations are in fact the only symmetry of this ring of monkey blocks --
it has none of the symmetries found in three dimensions.

\begin{wrapfigure}[13]{r}{0.3\textwidth}
\vspace{-10pt}
\centering
\includegraphics[width=0.3\textwidth]{Figures/monkeyblock_ring_left_screw}
\caption{A ring of four monkey blocks.}
\label{ring}
\end{wrapfigure}

Those familiar with four dimensional objects may notice that this ring of four cubes is half of the boundary of a hypercube. It might seem unlikely that by continuing to match the same monkey blocks face to face one could complete the hypercube without violating the matching rules. However,
not only does this work, but
the set of eight matched monkey blocks have symmetry around all three possible sets of isoclinic rotations: the $wx$-plane with the $yz$-plane, the $wy$-plane with the $xz$-plane, and the $wz$-plane with the $xy$-plane.
There are no other symmetries.

\begin{figure}[htbp]
\centering
\hspace{-.7cm}
\subfloat[Eight monkey blocks arranged into a hypercube net.]
{
\includegraphics[width=0.48\textwidth]{Figures/monkeyblock_hypercube_net_left_screw_right_mult_v3.jpg}
\label{monkeyblock_hypercube_net}
}
\hspace{.7cm}
\subfloat[The hypercube net labelled with the elements of $Q_8$. This can be interpreted as the Cayley graph for $Q_8$ with generating set $\{i,j,k\}$.]
{
\includegraphics[width=0.43\textwidth]{Figures/elements_left_screw_right_mult_v2}
\label{elements}
}
\caption{}
\label{monkey block hypercubes}
\end{figure}

In Figure \ref{monkeyblock_hypercube_net} the monkey blocks are shown as the net of a hypercube, with the three isoclinic rotations shown as pairs of axes. The monkey tails rotate around in two perpendicular loops of four blocks each, which we label as the $i$ axes. The $j$ axes are the two loops of four blocks made by connecting matching hands, and the $k$ axes go through the face faces. The $k$ axis is perhaps the easiest to see as the four blocks in one loop are arranged vertically, spinning up one step and a quarter turn to the left each time. At the same time, in the perpendicular plane the other four blocks move likewise: a quarter turn to the left around the vertical $k$ axis while rotating a quarter turn to the left around the horizontal $k$ axis.

Our labelling of the axes suggests that these isoclinic rotations can be thought of as multiplication by the quaternions $i$, $j$, and $k$, and of course this is true. Starting from an arbitrarily chosen block which we label 1, we can label each of the eight cells of the hypercube by the element of $Q_8$ we need to multiply (on the right) to get there. See Figure \ref{elements}.

Starting with any labelled cube, if we follow the arrows and write down the sequence of labels, then we finish at the cube labelled with the element of $Q_8$ which is the product of the starting cube label with the arrow labels, read in order. Following an arrow backwards corresponds to multiplying by the negative of the label on the arrow. Using these rules, we can see that rotating twice along any of the three axes (an isoclinic rotation of $180^\circ$) will take us to the opposite cube, as will rotating by $90^\circ$ once each along the three axes in order, geometrically showing that $i^2 = j^2 = k^2 = ijk = -1$.

\section{The quaternion group in sculpture}
\label{quat group in sculpture}

While the quaternion group is well known as an abstract group, as are its actions on objects such as the hypercube and the 3-sphere, to our knowledge no physical object has been made which has as its symmetry group the quaternion group.

It is a reasonably simple exercise to check that $Q_8$ is not a subgroup of the isometries of 3-dimensional Euclidean space. However, as noted in Section \ref{Sec:Quat group}, $Q_8$ is very naturally a subgroup of the symmetries of the 3-dimensional sphere, $S^3$, the unit sphere in $\mathbb{R}^4$. Using techniques we described previously in \cite{sculptures_in_S3}, we can create a sculpture that has $Q_8$ as its symmetry group. As noted in Section \ref{Sec:Quat group}, the eight elements of $Q_8$ act on the eight cubical cells of the hypercube. We radially project the boundary of the hypercube onto $S^3$, then use stereographic projection to move the eight cubes into $\mathbb{R}^3$, where we can build a physical sculpture using 3D printing. We must choose a design to put in each of the eight cubes. The design must itself have no symmetry, or else the symmetry group would be larger than $Q_8$. In order to be produced as a physical sculpture, the design in each cubical cell of the hypercube also needs to physically connect with its copies in the six neighbouring cubical cells, through the six square faces of the cube. The design must connect onto copies of itself following the same matching rules as the monkey blocks of Section \ref{Sec:Visualise quat group}.

\begin{wrapfigure}[26]{r}{0.63\textwidth}
\vspace{-5pt}
\centering

\includegraphics[width=0.63\textwidth]{Figures/euclidean_monkeys.jpg}

\caption{The monkey in a Euclidean cube, prior to mapping into $S^3$, with six neighbours.}
\label{euclidean_monkeys_in_cube}
\end{wrapfigure}

The obvious (perhaps even canonical) choice for such a design is, of course, a monkey. See Figure \ref{euclidean_monkeys_in_cube}. With appropriate posing of the monkey, its bilateral symmetry can be broken, and including the head and tail it has six limbs, one for each face of the cube. The monkey's left foot stands on the head of a neighbour, the left hand grabs a neighbour's right foot, and the right hand grabs a neighbour's tail. By symmetry, everything that goes around comes around -- so the other three neighbours of this monkey are standing on this monkey's head, grabbing its right foot, and grabbing its tail.

The monkey was designed in a Euclidean cube. It was then run through eight different transformations in order to move eight copies of it to the appropriate positions in $S^3$ and then back to $\mathbb{R}^3$ by stereographic projection. The first step of all of these transformations is to project the Euclidean cube into a curved cube in $S^3$. This is done in exactly the same way that a cubical cell of the hypercube is radially projected onto the hypersphere $S^3$. To be precise, we think of a point $(x,y,z)$ in the Euclidean cube $[-1,1]^3$ as actually being the point $(1,x,y,z)$ on one of the cells of the Euclidean hypercube $[-1,1]^4$, and map it to $S^3$ by

$$  (x,y,z)\mapsto \frac{(1,x,y,z)}{\sqrt{1+x^2+y^2+z^2}}.$$

Now that the design is on $S^3$, we (right) multiply it by $1,i,j,k,-1,-i,-j$ and $-k$ respectively for the eight transformations, and stereographically project each back to $\mathbb{R}^3$. There is a choice of where to project from --  we put the north pole at a vertex of the hypercube, so that in the projection the copies of the monkey are as far from infinity as possible. This makes the resulting features of the eight monkeys as near in size to each other as possible. Very small features may be too fragile to 3D print -- to avoid this we can scale the entire sculpture up, but only so far as the largest features fit within the printer and our budget.

The resulting sculpture is shown in Figure \ref{fig_monkeys}.  Note that the sculpture has no ``ordinary'' symmetries at all: every monkey is different if we only consider isometries of 3-dimensional space. However, under the appropriate isometries of the 3-sphere (as seen through the lens of stereographic projection) they are all identical. From this vantage point the three pairs of axes of rotation have equal billing: each circle consists of two larger, outer monkeys and two smaller, inner monkeys. The three pairs of axes go through the head-foot, hand-foot and hand-tail connections.

\begin{figure}[hp]
\centering

\subfloat[]
{
\includegraphics[width=.72\textwidth]{Figures/hypercube_monkeys_1_left_screw.jpg}
\label{monkeys_1}
}

\subfloat[]
{
\includegraphics[width=.72\textwidth]{Figures/hypercube_monkeys_2_left_screw.jpg}
\label{monkeys_3}
}
\caption{Two views of the sculpture ``More fun than a hypercube of monkeys''.
}
\label{fig_monkeys}
\end{figure}

Each monkey sits inside of a cell of the hypercube and connects to its neighbours through the 2-dimensional faces of the hypercube. Therefore, taken together they form the edges and vertices of the dual polytope to the hypercube, namely the 16-cell. Figure \ref{sixteencell} shows a 3D print of the edges and vertices of the 16-cell, projected to 3-dimensional space in the same way that the monkeys are. This perhaps shows the pairs of axes of rotation more clearly. The symmetry group of this object contains $Q_8$ but is much larger -- the asymmetric monkey design is necessary in order for us to represent \emph{only} the quaternion group.

\begin{wrapfigure}[10]{r}{0.3\textwidth}
\vspace{-10pt}
\centering

\includegraphics[width=0.3\textwidth]{Figures/updated_16-cell.jpg}

\caption{The 16-cell.}
\label{sixteencell}
\end{wrapfigure}

The object shown in Figure \ref{Fig:*44} is \emph{achiral}, meaning that it is the same as its mirror image. The object in Figure \ref{Fig:224} is \emph{chiral} -- the mirrored object is different.\footnote{To be precise, an object is \emph{achiral} if there is an orientation preserving isometry which moves the mirrored object to coincide with the original. An object is \emph{chiral} if this is only possible with an orientation reversing isometry.} However, the symmetry group of the mirrored object is the same as that of the unmirrored one.\footnote{Or rather, it is a conjugate of the original group of isometries by an orientation preserving isometry.} This is true of almost all objects, if we restrict ourselves to only consider the ordinary symmetries of 3-dimensional objects.\footnote{There are a few metachiral space groups in three dimensions, but no point groups.} However, our monkey sculpture with $Q_8$ symmetry is different. In Conway's classification it is \emph{metachiral} (see \cite{on_quaternions_and_octonions}, p44). Not only is it different from its mirror image, but also all of the left screw symmetries become right screw symmetries in the mirror image, and so the symmetry group itself is also reflected!\footnote{By which we really mean that it is a conjugate of the original group of isometries only by an orientation \emph{reversing} isometry.} In other words, the symmetry group itself has a handedness.

\section{Future directions}
\label{future directions}

It would be an interesting project to list the abstract groups against works of art or illustration, or objects from nature with the given group as its symmetry group.
Are there reasonably simple polytopes (perhaps shaped somewhat less like monkeys) that have symmetry group $Q_8$? Further new sculptures could be produced via stereographic projection, using other finite subgroups of the isometries of $S^3$, as listed in Tables 26.1 through 26.3 of \cite{symmetries_of_things}.
 Are there other natural strategies for making objects exhibiting previously unrepresented symmetry groups?

\section{Acknowledgements}

We thank Marc ten Bosch for motivating this line of thought with some interesting questions about computer graphics in four dimensions, and Will Segerman for design of the monkey.

\bibliographystyle{hyperplain}