\section{Introduction}
\noindent
In recent years, quantum annealers have received significant attention because of their potentials in solving complex computational problems which are often intractable with non-quantum computational methods~\cite{Farhi2001,Das2008,Albash2018,Hauke2020}. Quantum annealing is a procedure of making Hamiltonian $H(t)$ of a quantum many-body system adiabatically evolve from $\hat H_i$ to $\hat H_f$,
\begin{equation}
\hat{H}(t)=\hat{H}_i -\frac{t}{t_f}(\hat{H}_i-\hat{H}_f), \label{H1}
\end{equation}
so that the quantum state $\ket{\Psi(t)}$ initially prepared in the ground state of the former reaches the ground state of the later. Quantum annealers are considered with superconducting qubits~\cite{Johnson2011,Bunyk2014,Rosenberg2017,HarriganNP2021} and trapped-ion qubits~\cite{Korenblit2012,Islam2013,Richerme2013,Hauke2015}, aiming for various combinatorial optimization problems such as quantum simulations~\cite{Babbush2014}, classfications~\cite{Neven2008}, planning~\cite{Rieffel2015}, etc.  While many efforts in quantum annealing are being focused on large-scale implementations~\cite{Young2013,Boixo2014,Pudenz2014,Ebadi2020,Scholl2020} towards quantum speedup~\cite{Altshuler2010,Somma2012,Ronnow2014,Muthukrishnan2019}, here we explore the possibility of high-dimensional qubit connectivities. There are theoretical proposals emphasizing and thus utilizing qubit connectivities for NP-hard optimization problems: for examples, Lechner-Hauke-Zoller scheme~\cite{LHZ2015,Glaetzle2017} proposed a quantum annealing architecture for all-to-all connectivities using local interactions; and quantum optimization protocols are considered, e.g., for maximum independent set problems utilizing the nature of long-range couplings especially in Rydberg-atom quantum simulators~\cite{Pichler2018}. In the context relevant to the present work, Rydberg-atom quantum simulators~\cite{Weimer2010,Browaeys2020} draw attention, because of their high tunability in qubit connectivities~\cite{Lee2016,Kim2016,Barredo2016,Endres2016,Barredo2018} as well as many-body controllability in adiabatic processes~\cite{Pohl2010,Schauss2015,Bernien2017,Lienhard2018,Beterov2020}.

In this work, we consider, as a prototypical fractal structures, Cayley tree graphs of neutral atom arrangements, in which atoms and strongly interacting atom pairs respectively represent vertices and edges of the graphs (see Fig.~\ref{fig1}). Cayley trees are homogeneous and isotropic tree graphs of a fixed number of edges and no loop~\cite{Bethe1935,Ostilli2012}. Their infinite version is Bethe lattice widely used in various physics areas as a fundamental theoretical platform, often providing exactly solvable models in classical and quantum problems~\cite{Baxter1982}. In experiments to be described below, we construct atomic Cayley graphs of coordination number $Z=3$ and shell number up to $S=4$, as in Fig.~\ref{fig1}(a), and use a Rydberg-atom quantum simulator to directly probe their Ising-Hamiltonian phases.
\begin{figure*}[t]
\centerline{\includegraphics[width=1.00\textwidth]{Figure1.jpg}}
\caption
(a) A generic Cayley-tree graph of coordination number $Z=3$, with vertices and edges representing atoms and Rydberg-blockaded atom pairs, respectively. (b) A three-shell Cayley tree ($G_{10}$) of 10 atoms of inter-atom distance $d$. (c) A four-shell Cayley tree ($G_{22}$) constructed in three planes at $z=0, \pm h$ ($h=d/1.2$), where the last-shell branches are rotated by $72^\circ$ to avoid unwanted couplings. (d) A dual-center Cayley tree ($G_{14}$) of 14 atoms. (e,f,g) Plane-by-plane fluorescence images of the corresponding atom arrangements.}
\label{fig1}
\end{figure*}

\section{Cayley-tree atom arrangements in three-dimensional space}
\noindent
Neutral atoms (rubidium, $^{87}$Rb) are arranged in three-dimensional (3D) space with optical tweezers (far-off resonant optical dipole traps)~\cite{Barredo2018,Kim2020} (see Sec.~\ref{technical} for technical details). Three different Cayley-tree graphs are constructed. The first one is the three-shell ($S=3$) Cayley tree, which can be denoted by $G_{10}=(0s)(1s)^3(2s)^6$, having one in the first (center) shell, three in the second shell, and six in the third shell, as in Fig.~\ref{fig1}(b). The second one is the four-shell ($S=4$) Cayley tree, $G_{22}=(0s)(1s)^3(2s)^6(3s)^{12}$, having ten atoms in $G_{10}$ and twelve in the fourth shell, as in Fig.~\ref{fig1}(c).  And the last one is $G_{14}=(0s)^2(1s)^4(2s)^8$, a dual-center Cayley tree, having two first-shell atoms, four second-shell atoms, and eight third-shell atoms, as in Fig.~\ref{fig1}(d).

In our atom arrangements, each atom represents a vertex of a Cayley-tree graph and each pair of strongly-interacting atoms an edge. Non-connected atoms are supposed to interact with each other much weakly than connected atoms. The Hamiltonian of these atoms being coherently and simultaneously excited to a Rydberg energy state is given (in unit of $\hbar=1$) by
\begin{equation}
\hat{H}= \frac{1}{2} \sum_{j=1}^N\left\{\Omega\hat\sigma^{(j)}_x - \Delta\hat\sigma^{(j)}_z\right\}+ \sum_{j<k} U_{jk} \hat n^{(j)}\hat n^{(k)},
\label{H}
\end{equation}
where $N$ is the number of atoms, $\Omega$ is the Rabi frequency, $\Delta$ is the detuning, and $U_{jk}=C_6/|\vec{r}_j-\vec{r}_k|^6$ is the  pairwise atom interaction in the van der Waals interaction regime~\cite{Bernien2017,Lienhard2018}. Pauli operators $\hat\sigma_{x,z}$ are defined for a pseudo spin $1/2$ system composed of the ground state $\ket{\downarrow}=\ket{5S_{1/2}, F=2, m_f=2}$ and Rydberg state $\ket{\uparrow}=\ket{71S_{1/2},F'=3,m_F'=3}$ of each atom, and $\hat{n}=(1+\hat\sigma_z)/2$. For the strong and weak interactions of connected and nonconnected atom pairs, respectively, we set the distances of connected atoms the same and within the Rydberg-blockade radius, i.e., $|\vec{r}_j-\vec{r}_k|=d<r_b\equiv (C_6/\hbar\Omega)^{1/6}=9.8~{\rm \mu m}$, for $(j,k)\in E$ (the edge set of a graph $G$), and the distances of all others are $|\vec{r}_j-\vec{r}_k|>r_b$.

The two-dimensional (2D) arrangement of $G_{10}$ is made with $d<r_b<\sqrt{3}d$ to satisfy the above condition, in which, as shown in Fig.~\ref{fig1}(e), the last-shell atoms of different branches are separated more than $\sqrt{3}d$, having at most 1/27 times smaller  interactions than $U\equiv C_6/d^6$ of the connected atom pairs.
However, $G_{22}$ and $G_{14}$ cannot be planar, because in 2D arrangements, atoms of different outermost branches are too close, requiring nonplanar, three-dimensional (3D) arrangements. As shown in Fig.~\ref{fig1}(c), we rotate the last-shell branches of $G_{22}$ by an angle of $2\pi/5$ about the axes along the previous branches, so that all the last-shell atoms are well separated (more than $\sqrt{3}d$) from each other. In $G_{14}$, some last-shell branches are rotated similarly, as shown in Fig.~\ref{fig1}(d). As-constructed atom arrangements are shown in Figs.~1(e-g) respectively for $G_{10}$, $G_{22}$, and $G_{14}$.

\section{Phase diagrams of Cayley-tree Ising spins}
\noindent
With atoms arranged on one of the above Cayley-tree graphs, we perform quantum annealing to find the ground state of a target Hamiltonian. The atoms are initially prepared in $\ket{\downarrow\downarrow\cdots\downarrow}$ (the paramagnetic down spins in Phase I) under the Hamiltonian conditions of  $\Delta \gg U$ and $\Omega=0$ in $\hat H$. The target (final) Hamiltonian  is
\begin{eqnarray}
\hat{H}_G (U,\Delta_f)=U \sum_{(j,k)\in E}\hat n^{(j)}\hat n^{(k)} -\frac{\Delta_f}{2} \sum_{j=1}^N \hat\sigma^{(j)}_z
\label{HG}
\end{eqnarray}
where $E$ is the edge set of $G\in \{ G_{10}, G_{22}, G_{14}\}$, $\Delta_f$ is the final detuning, and no couplings are assumed for unedged atom pairs. We note that $\hat H_G$ is an Ising spin-glass Hamiltonian, given by
\begin{eqnarray}
\hat{H}_G = J \sum_{(j,k)\in E} \hat \sigma_z^{(j)}\hat \sigma_z^{(k)}
+   h_z^C    \sum_{j\in C}  \hat \sigma^{(j)}_z + h_z^V   \sum_{j\in V }  \hat \sigma^{(j)}_z,
\label{HG2}
\end{eqnarray}
where $J=U/4$ is the coupling, $h_z^C= 3U/4-\Delta_f/2$ and $h_z^V=  U/4-\Delta_f/2$ are the local fields, and $C$ and $V$ denote the vetices in the core (inner shells) and valence (outermost) shell, respectively, of a Cayley-tree graph $G(E,C+V)$. In general, as the ground-states of an Ising Hamiltonioan depend on the specific atom arrangement of $G$, finding such arrangement-specific ground-states of an arbitrary Ising-spin graph is an NP-complete problem~\cite{Lucas2014}.

However, Cayley-tree graphs allow heuristic understanding of their phase diagrams as follows: Cayley trees have more atoms on the valence shell than on the rest (the inner shells) of a tree. Therefore, with Hamiltonian $\hat H_G$ in Eq.~\eqref{HG} (of positive coupling, $U>0$, in our consideration), the valence spins are all aligned either up or down depending upon $\Delta_f>0$ or $\Delta_f<0$, respectively, resulting in the $\Delta_f=0$ phase boundary. For $\Delta_f<0$, inner-shell spins, adjacent to the valence spins, also favor down spins, as their couplings, $U\hat n^{(j)} \hat n^{(k)}$, to adjacent outer-shell down spins are zero, which sequentially results in all down spins (so paramagnetic down spins, Phase I). For $\Delta_f>0$, the single-spin flipping energy ($\Delta_f$) competes with adjacent anti-ferromagnetic couplings. For $\Delta_f>3U$, the former is always higher than the latter (of max $3U$), so all spins are up (paramagnetic up spins, Phase II). In between, $0<\Delta_f<3U$, the anti-ferromagnetic coupling is stronger, thus favoring the anti-ferromagnetic phase (Phase III) of shell-by-shell alternating spins. The resulting $\hat H_G$ phase diagram is shown for $G_{10}$ in Fig.~\ref{fig2}(a), in which the anti-ferromagnetic phase (Phase III) has up spins on the first and thrid shells and down spins on the second. The phase diagram for $G_{22}$ differs from $G_{10}$ only in Phase III, where the odd-number shells are down spins and the even-number shells are up spins. In short, the phase diagram of $\hat H_G$ (of positive $U$) for regular Cayley-graphs ($G_{10}$ and $G_{22}$ in our consideration) has the para-antiferro phase boundaries determined by $\Delta_f=3U$ and $\Delta_f=0$, so the anti-ferromagnetic ordering is expected in the region $0<\Delta_f<3U$ (Phase III), as in Fig.~\ref{fig2}(a).

\begin{figure}[tb]
\centerline{\includegraphics[width=0.50\textwidth]{Figure2}}
\caption{Phase diagrams of Cayley-tree Ising Hamiltonians (a) $\hat H_{G_{10}}$ and $\hat H_{G_{22}}$ and (b) $\hat H_{G_{14}}$, in which ground-state spin configurations are paramagnetic down spins (Phase I), paramagnetic up spins (Phase II), antiferromagnetic phase (Phase III), and antiferro-like phases with center-spins $\ket{\uparrow\uparrow}$ (Phase IV) and $(\ket{\uparrow\downarrow+\downarrow\uparrow})/\sqrt{2}$ (Phase V).
Circled numbers indicate the Hamiltonian parameters of the experimental data in Figs.~\ref{fig3} and \ref{fig4}: {\textcircled{\small 1}} ($U/\Omega_0$, $\Delta_f/\Omega_0$)=(1.82, 2) in Phase III of $G_{10}$,
{\textcircled{\small 2}} (2.25, 2) in Phase III of $G_{22}$,  {\textcircled{\small 3}} (1.67, 2) in Phase IV of $G_{14}$,
{\textcircled{\small 4}} (2.70,2) and {\textcircled{\small 5}} (5.41,2) in Phase V of $G_{14}$. $\Omega_0=1.1$~$(2\pi)$MHz in all experiments.
}
\label{fig2}
\end{figure}

The phase diagram for $G_{14}$, the dual-center Cayley tree, is a little more complex than regular Cayley trees, because of the frustrations of the center spins. Exact diagonalization of $\hat H_{G_{14}}$ finds the phase diagram as shown in Fig.~3(b), in which the paramagnetic phases (Phases I and II) are the same as $G_{10}$, but Phase III is split to Phases IV and V, of respective center-spin configurations of $\ket{\uparrow\uparrow}$ and $(\ket{\uparrow\downarrow+\downarrow\uparrow})/\sqrt{2}$. The new phase boundary is $U=\Delta_f$, along which the single-spin flipping cost  equals the additional frustration cost of the center spins. It is noted that, in Phase V, other energy-degenerate ground states, with respective center-spin configurations of $\ket{\uparrow\downarrow}$, $\ket{\downarrow\uparrow}$, and $(\ket{\uparrow\downarrow-\downarrow\uparrow})/\sqrt{2}$, are all dark states for a quantum annealing performed to the given initial state, due to the symmetey of $\hat H$ in Eq.~\eqref{H}.

\section{Experimental verification of Cayley-tree Ising Phases}
\noindent
Experiments are performed to verify the antiferro-like ground states (the Phases III, IV, and V) of Cayley-tree Ising spins (see Sec.~\ref{technical} for technical details). Quantum annealing is proceeded with three stages of time evolution, along a vertical control path from Phase I to either Phase III, IV or V in Fig.~\ref{fig2}. In the first stage ($0<t<0.1t_f$), Rabi frequency is adiabatically turned on from $\Omega_i=0$ to $\Omega_0=1.1 (2\pi)$MHz, while detuning is maintained at $\Delta=-2\Omega_0$ for the paramagnetic down-spin ordering ($|h_z^{C,V}| \gg J$) of Phase I. In the second stage ($0.1<t/t_f<0.9$), the detuning is swept from $-2\Omega_0$ to $2\Omega_0$, while the Rabi frequency is maintained at $\Omega_0$. In the final stage ($0.9<t/t_f<1$), the detuning is maintained at $\Delta_f=2\Omega_0$ (for anti-ferromagnetic ordering, $0<\Delta_f<3U$) and the Rabi frequency is adiabatically turned off from $\Omega_0$ to $\Omega_f=0$. The total operation time $t_f=(2\pi)3.2/\Omega_0=2.9$~$\mu$s is chosen long enough to identify the max-population states and short compared to the coherence time of 10~$\mu$s~\cite{Kim2020}. After the time evolutoin, the resulting spin configuration is detected and the procedure is repeated until the probability distribution of all spin configurations is obtained.
The final Hamiltonians $\hat H_G(U,\Delta_f=2\Omega_0)$ are chosen at various phases in Fig.~\ref{fig2}: {\textcircled{\small 1}} and {\textcircled{\small 2}} for Phase III of $G_{10}$ and $G_{22}$, respectively; {\textcircled{\small 3}} for Phase IV of $G_{14}$; and {\textcircled{\small 4}} and {\textcircled{\small 5}} for Phase V of $G_{14}$.
\begin{figure*}[t]
\centerline{\includegraphics[width=1.00\textwidth]{Figure3.jpg}}
\caption
Probability distribution of all spin configurations in enumerated bare-spin basis: (a) $G_{10}$ Cayley-tree atoms measured at {\textcircled{\small 1}} in Phase III and accumulated over 672 events. (b) $G_{22}$ at  {\textcircled{\small 2}} in Phase III (2208 events). (c) $G_{14}$ at  {\textcircled{\small 4}} in Phase V (5113 events). Max populations are ground states (highlighted in red) and smaller peaks are identified near-ground excited states. (d,e,f) Fluorescence images of ground-state ($\ket{\downarrow}$) atoms in characteristic spin configurations.}
\label{fig3}
\end{figure*}
\noindent

Measured probability distributions are shown in
Figs.~\ref{fig3}(a-c), which plot the probability distributions, respectively measured at {\textcircled{\small 1}}, {\textcircled{\small 2}}, and {\textcircled{\small 4}}, for all spin configurations of the Cayley-tree Ising spins. Spin configurations are represented in the bare-atom basis with binary enumeration of  $\ket{\downarrow}=\ket{0}$ and $\ket{\uparrow}=\ket{1}$~\cite{Kim2018}. After the quantum annealing, the wavefunction of the atoms is expectedly driven near the ground state of $\hat H_G(U, \Delta_f)$, so the corresponding spin configuration is measured with a high probability. In Fig.~\ref{fig3}(a), Phase III of $G_{10}$ is probed at {\textcircled{\small 1}}. The observed max-population state is $\ket{\uparrow;\downarrow\downarrow\downarrow;\uparrow\uparrow\uparrow\uparrow\uparrow\uparrow}=\ket{2^{9}+2^5+2^4+2^3+2^2+2^1+2^0}=\ket{575}$, which is the anti-ferromagnetic phase of shell-by-shell alternating spins, agreeing with the expected ground state of $\hat H_{G_{10}}$ in Phase III. Likewise, in Figs.~\ref{fig3}(b), Phase III of $G_{22}$ is probed at {\textcircled{\small 2}} and the observed $\ket{\downarrow;\uparrow^3;\downarrow^6;\uparrow^{12}}=\ket{1839103}$ agrees with the ground-state of $\hat H_{G_{22}}$ in Phase III. In Fig.~\ref{fig3}(c), which shows the measurement of $G_{14}$ at  {\textcircled{\small 4}} in Phase V, high populations are observed for $\ket{\downarrow\uparrow;\downarrow^4;\uparrow^8}=\ket{4351}$ and $\ket{\uparrow\downarrow;\downarrow^4;\uparrow^8}=\ket{8447}$, agreeing with the ground superposition state (of center spins in $(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow})/\sqrt{2}$) in Phase V of $\hat H_{G_{14}}$. Atom images (the fluorescence of ground-state atoms) of characteristic spin configurations are shown in Figs.~\ref{fig3}(d,e,f), including the max-population states (of star marks).

Phase IV of $G_{14}$ is also probed acrossing the IV-V phase boundary, with three different Cayley trees of respective edge lengths $d/r_b=0.76$, 0.86, and 0.92, which correspond to {\textcircled{\small 5}}, {\textcircled{\small 4}}, and {\textcircled{\small 3}} in the phase diagram. In Fig.~\ref{fig4}, the probabilities of three high-population states, of respective center-spin configurations $\ket{\uparrow\uparrow}$, $\ket{\downarrow\downarrow}$, and $(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow})/\sqrt{2}$, are plotted. The max-populated state changes from $\ket{\uparrow\uparrow}$ ($d<d_c=0.89r_b$, Phase IV) to $(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow})/\sqrt{2}$ ($d>d_c$, Phase V), agreeing with the phase boundary given by $U(d_c)=\Delta_f$. It is noted that the higher-order long-range couplings which are ignored in Eq.~\eqref{HG} play little role in the tested parameter region and that other ground states, of center spins $\ket{\downarrow\uparrow}$, $\ket{\uparrow\downarrow}$, and $(\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow})/\sqrt{2}$) are forbidden by the Hamiltonian symmetry. Also, the first excited state changes from $(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow})/\sqrt{2}$ to $\ket{\downarrow\downarrow}$ and then to $\ket{\uparrow\uparrow}$ from left to right, in accordance with numerical calculation.

\begin{figure}[tb]
\centerline{\includegraphics[width=0.50\textwidth]{Figure4}}
\caption{Probabilities of the low-energy states of Cayley-tree Hamiltonian, $\hat H_{G_{14}}(U(d), \Delta_f=2\Omega_0)$, measured for three different edge lengths, $d/r_b=0.92$, 0.86, and 0.76, which correspond to  {\textcircled{\small 3}}, {\textcircled{\small 4}}, and {\textcircled{\small 5}}, respectiveily, in the $G_{14}$ phase diagram in Fig.~\ref{fig2}(b).}
\label{fig4}
\end{figure}

\section{N\'{e}el's order formation dynamics}
\noindent
Furthermore, we measure the phase formation dynamics during the quantum annealing process. Figure~\ref{fig5} shows the time evolution of N\'{e}el's order, defined by $O_N\equiv -\sum_{(i,j)\in E} \langle\hat\sigma^{(i)}(t) \hat\sigma^{(j)}(t)\rangle /||E||$, where $||E||$ is the number of edges, along with the up-spin probabilities of individual atoms in $G_{10}$. With the snap-shot measurements, the adiabatic order formation is clearly observed from the initial paramagnetic phase to the final anti-ferromagnetic phase. The oscillatory behavior is attributed to the finite size. In comparison, a numerical calculation (solid lines) is performed with Lindblad master equations, taking into account state-preparation-and-detection (SPAM) errors ($P(\ket\downarrow|\ket\uparrow)=0.18$, $P(\ket\uparrow|\ket\downarrow)=0.02$), individual dephasing ($\sim$36~kHz $\ll \Omega_0$) due to the spontaneous decay during Rydberg excitations, and collective dephasing ($\sim$3~kHz) from laser phase noise~\cite{Lee2019}. The numerical calculation in accordance with the observed maximal N\'{e}el's order of $O_N (t_f) = 0.48(2)$ indicates that errors are largely due to the SPAM errors accumulated for $N=10$ atoms. The calculation also suggests that the maximal N\'{e}el's order before measurements was $O_N(t_f)=0.82$ and that the ground-state probability was 61\% through the quantum annealing.

\begin{figure}[t]
\centerline{\includegraphics[width=0.49\textwidth]{Figure5}}
\caption
Quantum annealing dynamics of the N\'{e}el's order in the $G_{10}$ Cayley tree are probed as a function of the evolution time and compared with numerical calculations. At chosen times, the Rydberg-state probabilities, $\langle{\hat n}\rangle_j$, of all atoms ($j =1, \cdots, 10$) are shown, being plotted at their respective atom sites.}
\label{fig5}
\end{figure}

\section{Experimental details} \label{technical}
\noindent
The above experiments were performed with a Rydberg-atom quantum simulator previously reported elsewhere~\cite{Kim2018,Lee2019,Jo2020,Kim2020}. In the quantum simulator, rubidium ($^{87}$Rb) atoms were initially prepared in the hyperfine ground state, $\ket{\downarrow}=\ket{5S_{1/2},F=2, m_F=2}$, and trapped with optical tweezers. A spatial light modulator (SLM, Meadowlark ODPDM512) was used to create the 3D array of $2N$ optical tweezers and an electrically focus-tunable lens (EL-16-40-TC from Optotune) verified the positions of captured atoms in each atom plane. The SLM was computer-programmed with weighted Gerchberg-Saxton (w-GS) algorithm~\cite{DiLeonardo2007,Kim2019} so that the resulting electric field near each target site $(x, y, z)$ was created as $E(x, y, z) = E_0 \sum_{X,Y}e^{i\Phi(X,Y)}e^{-iT}$, where $X,Y$ are the SLM coordinates on the Fourier plane of the optical tweezers, $\Phi(X,Y)$ is the SLM phase pattern, and $T= {2\pi (x X+y Y)}/{f\lambda } + {\pi z(X^2+Y^2)}/{ f^2 \lambda}$ is the transfer kernel. The phase pattern $\Phi$ (for $2N$ optical tweezers) was obtained with the w-GS condition, given by $\Phi(X,Y)={\rm arg}(\sum_{j=1}^{2N} {w_j E_j(x,y,z)}e^{iT_j}/ {\abs{E_j(x,y,z)}})$, where each weighting factor $w_j$ for $j^{\rm th}$ optical tweezer was optimized through adaptive iterations. The number of the w-GS iterations was about 5, taking about 20~ms time to generate an array of 50 tweezers. Determinstic rearrangement of $N$ atoms to target sites, about 20~$\mu$m apart from the initial atom reservoir, was performed with consecutive 45 frames of moving traps programmed with phase induction~\cite{Kim2019}. Over 90\% target-occupation probabilities were achieved in 900-ms reconfiguration time and all measurements were performed with defect-free arrangements.

After the atom array was prepared, the optical tweezers were temporarily turned off and the quantum annealing was proceeded. The atoms were excited to the Rydberg state, $\ket{\uparrow}=\ket{71S_{1/2},m_J=1/2}$, via the off-resonant intermediate state, $\ket{i}=\ket{5P_{3/2}, F'=3, m_F'=3}$. Two lasers (780-nm and 480-nm lasers for $\ket{\downarrow}\rightarrow\ket{i}$ and $\ket{i}\rightarrow\ket{\uparrow}$) were used for the two-photon transition~\cite{Kim2020}.  $\Omega(t)$ and $\Delta(t)$ in $H(t)$ were programmed with a RF synthesizer (Moglabs XRF, 10~MHz) of frequency (780 nm) and amplitude (480 nm) modulations. The modulation ranges were $0\le\Omega(t)\le\Omega_0$ and $-2\Omega_0\le\Delta(t)\le2\Omega_0$, where the max Rabi frequency was $\Omega_0=\Omega_{780}\Omega_{480}/(2\Delta')=1.1$~$(2\pi)$MHz given by $\Omega_{780}=75$~$(2\pi)$MHz, $\Omega_{480}=19$~$(2\pi)$MHz, and $\Delta'=660$~$(2\pi)$MHz (the intermediate detuning). The van der Waals coefficient~\cite{Weber2017} was $C_6=(2\pi)1004$~GHz$\times{\mu m^6}$ for $\ket{71S}$ Rydberg-state atoms and the Rydberg blockade radius was $r_b=({C_6/\Omega_0})^{1/6}=9.8~{\rm \mu m}$.
The frequency error due to AC Stark shift was small, below 140 kHz, mainly caused by the 480-nm laser.  After the quantum annealing, the optical tweezers were turned back on and the atom states in the ground state ($\ket \downarrow $) were measured whether they survived ($\ket \downarrow $) or not ($\ket \uparrow$). Each measurement in Fig.~\ref{fig3} was repeated by 672, 2208, and 5113 times respectively for $G_{10}$, $G_{22}$, and $G_{14}$ to obtain the accumulated probability distributions.

\section{Conclusion}
\noindent
We have explored the possibilites of Rydberg-atom quantum simulators towards high-dimensional qubit connection programming. With up to $N=22$ rubidium single atoms arranged in three-dimensional space, we have programmed Ising Hamiltonians on three different Cayley tree graphs. The anti-ferromagnetic phase in regular Cayley trees and frustrated competing ground-states in a dual-center Cayley tree are directly observed, showing good agreement with model calculations. There are many hopes to utilize the potentials of Rydberg-atom quantum simulators in large-scale implementations and high-dimensional configurability~\cite{LHZ2015, DiehlNP2011, Glaetzle2017, Pichler2018, Serett2020,Qiu2020}. It is hoped that 3D-qubit configurations of Rydberg-atom quantum annealers
shall be useful for efficient and programmable quantum optimization problems.

\begin{acknowledgements}
This research was supported by Samsung Science and Technology Foundation (SSTF-BA1301-52), National Research Foundation of Korea (NRF) (2017R1E1A1A01074307), and Institute for Information \& Communications Technology Promotion (IITP-2018-2018-0-01402).
\end{acknowledgements}


\section{Introduction}
\label{sec1}
Colourisation is the process of inferring colour information from a grey-scale prior and then combining the colour, with the prior, to produce a full-colour image. It is an ill-posed problem. For every prior, there are many plausible colourisations and many more that are implausible. Manual colourisation, as carried out by humans, relies on heuristics, common sense and prior knowledge. There is a large body of work on semi-automatic and automatic colourisation. The most recent work in colourisation involves the use of Generative Adversarial Networks (GANs) \cite{goodfellowGAN2014}, in which one deep neural network ($\mathcal{G}$) tries to create plausible colourisations and a second network ($\mathcal{D}$) tries to discriminate between real colourisations and those produced by $\mathcal{G}$. For the colourisation task, GANs try to learn the manifold on which all the plausible colourisations lie. For training, a set of natural images is converted to a colour space that separates the grey-scale from the colour information; CIEL*a*b* is the most popular for this. $\mathcal{G}$ is then trained to produce a*b* from an L* prior. This training regime, only admits a single plausible colourisation for each prior, which is limiting, when trying to learn the full manifold of all plausible colourisations. There have been some attempts to measure the performance of colourisations, but all are problematic, particularly when it comes to applying a similar score to multiple plausible colourisations. This paper will present one possible measurement system.

\section{Measuring Colour}
Most colourisation techniques rely on some form of human visual inspection to determine efficacy, or for comparison to other techniques, \cite{IizukaLetThereBeColor:2016:LCJ:2897824.2925974}, \cite{IronyColorByExample2005}, \cite{Isola2017}, \cite{WelshTransferingColor:2002:TCG:566654.566576},  \cite{YatzivChrominanceBlending1621234},  \cite{ZhangIE16ColorfulImageColorization}, \cite{ZhangZIGLYE17UserGuidedColorization}.  Human visual inspection can include qualitative analysis, naturalness scoring, or real/fake judged by human. In some cases the Amazon Mechanical Turk (AMT) was used with up to 40 participants looking at each technique, \cite{ZhangIE16ColorfulImageColorization}.

Many attempt an objective measure based on pixel value errors, such as RMS pixel error, \cite{Deshpande7410429LargeScaleAutoColorization}, PSNR \cite{Cheng2015DeepC},  \cite{ZhangZIGLYE17UserGuidedColorization}, L2 pixel distance  \cite{ZhangIE16ColorfulImageColorization} (only within a defined threshold).  Some methods \cite{ZhangIE16ColorfulImageColorization}, \cite{larsson2016learning}, \cite{Vitoria_2020_WACV}  utilise the concept that colour assists in classifying objects. Therefore a network designed to classify objects using colour images will show a deterioration in performance if inferred with a poorly colourised image. The difference can then be used as a proxy measure for colourisation performance.
From the results of \cite{ZhangIE16ColorfulImageColorization}, it is clear that while correct colour images achieve higher classification scores than their equivalent grey-scale image, the grey-scale image still makes up the majority of the information required for classification. E.g. grey-scale 52.7\% vs ground-truth colour 68.3\%.

As GANs are now a primary focus of colourisation research, it is worth mentioning two measures, Fr\'echet Inception Distance, \cite{HeuselRUNKH17} and Inception Score, \cite{SalimansNIPS2016_6125}, that are starting to gain a consensus in GAN literature. However, it is clear that they are imperfect measures and have not, to our knowledge, been used for the judgement of the fidelity of colourisations.

\subsection{Problems in measuring colourisations}
The first problem with measuring colourisation is that it is a multi-modal problem. While those colourisation methods that employ hints, from the user, require a single colour target for each grey-scale prior, fully-automatic colourisation allows for many plausible colourisation targets.

The second problem is that of spatial coherency. The colour should remain consistent within a given semantic region, and then change precisely at the boundary between regions.

Pixel difference methods (L1, L2, PSNR, MSE, RMS), measure only the difference at a pixel, in isolation from its neighbours, between a colourisation and a single ground-truth. So they penalise any sample from the distribution of plausible colourisations that is not the ground-truth . They fail to penalise for poor spatial coherency as they don't take account of neighbouring pixels, \cite{MulleryWhelan2019}.
This has led to the wide-spread use of naturalness measures based on human inspection.

A third problem, which we do not consider in this work, is the plausibility of the modes of colour for a given semantic object, e.g. Blue and Green are implausible colours for skin tones, but plausible for clothing.

\subsection{What makes a plausible colourisation?}
As stated above, we will not look at the plausibility of colour modes for specific semantic objects. The section on future work will outline potential routes forward which may benefit from the work in this paper. Instead, we will make the simplifying assumption that plausible colourisations can take on any hue, as long as that hue is consistent within the borders of the semantic object.
We define consistent as maintaining the same small value changes (texture) in hue as the ground-truth within the borders of a semantic segment and changing the hue value by a significant amount at the same spatial locations as these changes in a ground-truth image.

The saturation/chrominance of an object depends on the object itself, its reflectance properties (specular/lambertian) and the lighting conditions. Human colourists consider all of these things when deciding what level of saturation to apply. With this in mind, we can see that each object has a narrower distribution of plausible saturations that could be applied, compared to hue. We make the simplifying assumption that the saturation of an object should be the same in all colourisations as in the ground-truth.
It is difficult to make such simplifying assumptions about the CIEL*a*b* space, as equal chrominance (a close relative of saturation) forms a circular path around the centre of a*b* space.
However, CIEL*a*b* attempts to provide a coordinate system between two stimuli that is close to the colour difference perceived by the human eye over small distances, \cite{BRAINARD2003191}, and therefore we will convert CIEL*a*b* to $L^*h_{ab}^{*}c_{ab}^{*}$, \cite{fairchild.ch10}, using the following equations.
\begin{equation}
c_{a b}^{*}=\sqrt{\left(a^{* 2}+b^{*} 2\right)}
\end{equation}
\begin{equation}
h_{a b}^{*}=\tan ^{-1}\left(b^{*} / a^{*}\right)
\label{hab}
\end{equation}
For the rest of this paper, $h_{a b}^{*}$ will be referred to as hue, and $c_{a b}^{*}$ will be referred to as chroma.

\section{Proposed method}
We propose to measure the chroma using the Structural Similarity Index Measurement (SSIM), \cite{SSIM}. SSIM is designed to measure the perceived quality. It takes into account local spatial comparisons of magnitude/luminosity, contrast and structure. We give each of these equal precedence, as did the original authors of SSIM (see equation (\ref{SSIMComb}) below, $\alpha=\beta=\gamma=1$).
\begin{equation}
SSIM(x, y)=[l(x, y)]^{\alpha} \cdot[c(x, y)]^{\beta} \cdot[s(x, y)]^{\gamma}
\label{SSIMComb}
\end{equation}

We cannot use SSIM to measure hue for the following reasons. Firstly, the magnitude/luminosity measure of SSIM will penalise in the same way as the pixel difference measures, and this means those plausible hues that are not close to the ground-truth will be penalised.
Secondly, the local structure measurement is also problematic in that it requires the sign of local changes to be identical between the images being compared. Local changes in hue, for two plausible colourisations, would not need to be of the same sign.

The contrast measurement of SSIM is more promising, as two plausible colourisations might be expected to have the same contrast in hue within a semantic segment of the image. A problem occurs at the border of two segments. The contrast here will be relatively large, but the actual magnitude of the difference would not necessarily be the same between two plausible colourisations. So our third issue, with SSIM for hue, is that it cannot treat all large changes in hue as being similar, which is necessary to compare plausible colourisations.

\subsection{Other nuisances of hue space}
Hue represents an angle. There is no hard edge in the CIEL*a*b* colour space. Hue here should not be mistaken with the idea of a specific frequency of colour. Colour frequencies do not rotate back circularly. However, from the perspective of human perception, a combination of Red and Blue (opposite ends of the spectrum), is possible and as such the hue does complete a continuous periodic space. If two hues are quite close together and so may represent a colour texture, they may be represented by numbers that are far apart. Using equation (\ref{hab}) we will obtain hue values $\in [-\pi,\pi]$. However, in our space, $-\pi$ and $\pi$ are zero distance apart. For any number system we choose, these will be considered the maximum distance apart. We instead, ensure that the shorter distance around the circle is always chosen with what we will refer to as an \textit{Angle\_Diff} function, see equation (\ref{angleDiff}).

A second problem to consider is the hue of negligible chroma values. In circumstances where the chroma value is minimal, then the specific hue has little meaning. Even large changes in hue will not be noticeable to a viewer if the associated chroma at that location is very small. Despite this, we have no mechanism to convey an absence of hue. In many cases, the hue can be noisy and have changes that appear significant when analysed without consideration of the chroma. This occurs when negligible values of a* and b* can be modified by small rounding errors in calculation, often leading to large changes in hue using equation (\ref{hab}). There is no added utility for a colourisation algorithm to attempt to generate these invisible changes in hue. In our proposed algorithm, we will multiply the gradients of the hue channel by a logistic function of the chroma value for each pixel location to disregard gradients in the hue that corresponds to negligible chroma values.

\subsection{Spatial Hue Similarity Measure SHSM}
\begin{equation}
	\nabla h'= \mbox{\textit{Angle\_Diff}}\left(\left|\nabla  h\left(\frac{1}{1+e^{-k\left(c-c_{0}\right)}}\right)\right|\right)
	\label{angleDiff}
\end{equation}

$h$ is the hue channel, and $c$ is the chroma channel. $c_0$ is the user selected centre of the chroma range that we use for the logistic function to reduce sensitivity to hues that correspond to minimal chroma values. $k$ is the slope of the curve at its mid-point.
As hue is a circular space, the \textit{Angle\_Diff} function is used to find the smallest gradient that represents the change in hue. We use the absolute value of hue change $(| |)$ to be insensitive to the direction of change in the hue space.

\begin{equation}
\nabla h_n=\frac{1}{1+e^{-k\left(\nabla h'-h_{0}\right)}}
\label{logistic}
\end{equation}

Here we transform the $\nabla h'$ through another logistic function. $h_{0}$ is the chosen mid-point of the curve and $k$ is the slope of the curve at its mid-point. The output of this function will be $\mathbb{R} \in [0,1]$. We denote this, $\nabla h_{n}$ as it is the hue gradients normalised to the range $[0,1]$.

Transformation, through the logistic function, attenuates the significance of the smallest hue changes that are less perceptible to the human eye. It also extends the significance of small to medium changes in the hue that would represent perceptually visible hue texture. Finally it constrains, to a value of $1$, all hue changes in the medium-to-large range that would be consistent with a change in the semantic region.

Just as in SSIM, to determine spatial similarity we will convolve ($\circledast$) $\nabla h_{n}$ with a 2D Gaussian $g$, we retain the values for this Gaussian as used in the original SSIM default formulation, i.e. $\sigma=1.5$, \cite{SSIM}.

\begin{equation}
	\nabla h_{s} = g \circledast \nabla h_n
	\label{}
\end{equation}
Where $\nabla h_{s}$ represents the gradient of spatially weighted hue.

Finally we use the contrast part of SSIM to do our hue comparison.
\begin{equation}
	\textit{SHSM}(GT, CL) = \frac{2 \nabla h_{s}(GT)\nabla h_{s}(CL) + \epsilon}{\nabla h_{s}(GT)+\nabla h_{s}(CL) + \epsilon}
	\label{hueContEq}
\end{equation}

Here, $GT$ is the ground-truth image, and $CL$ is the colourisation that is being compared. As in SSIM, $\epsilon$\footnote{In the original SSIM paper this was called C, but we do not want this to confuse our notation where c is chroma.} is included for numerical stability but can also be used to softly ignore small, less perceptible, differences between the two images. We choose a value of 0.03\footnote{In SSIM, this was the value for K and C was derived from that with C=L*K where L was 255 for an 8-bit image. However, as we have normalised all values to 1, then C and K are the same value}, which is what SSIM chooses for the contrast part of the calculation.
We must calculate the gradients in two orthogonal directions (we choose the two diagonals), and so we will calculate $\textrm{SHSM}$ for the two directions and multiply the results.

\subsection{Hyper-parameter Choices}
To analyse the hue channel, 4000 natural images from the Places dataset, \cite{Places}, are converted to the $L^*h_{ab}^{*}c_{ab}^{*}$ space. The gradient of the hue channel is what is of interest, so a histogram of gradients in two orthogonal directions is constructed.
We normalise the histogram to show the percentage of pixels on the y-axis, see Fig \ref{removeNuicance}. The graph labelled 'Unprocessed' has many large spikes throughout the range, showing that specific gradients are very prevalent throughout the dataset. Each of these can be explained by a specific case of changes in small values of chroma. e.g. when converting from CIEL*a*b*, if $(a=0,b=0)$ changes to $(a=0,b=1)$ or $(a=1,b=0)$ we get a hue of $|\pm\frac{\pi}{2}|$ which is 63 when represented with 8-bits covering $-\pi \textrm{ to} +\pi$. This case represents the peak visible in the centre or the distribution. Similar peaks occur for combinations of $a$ and $b$ for values $0-6$. These do not represent accurate hue values but are instead a nuisance of the conversion to a polar form and as the chroma is a small value the hue change will be difficult to perceive. To remove the influence of these, we pass the hue gradients through a logistic function, Equation (\ref{angleDiff}). To choose an appropriate value of $c_0$ we reprocess the dataset with values of $c_0 \in [1,10]$ and re-construct the histogram, see Fig \ref{removeNuicance}.
 Increasing values of $c_0$ smooth the histogram.

We choose a value of $c_0=5$, as this appears to be the first of the lines to lose peaks in the distinctive locations of small chromas.

With these nuisances removed using a fixed value of $c_0=5$, we re-analyse the 4000 images to see how much variance is in the individual gradient histograms of the set, see Fig \ref{textureVBoundary}.
This provides valuable insight into the statistics of the hue channel. The majority of pixels have a gradient of $\leq2$. Gradients of this size will be challenging for humans to perceive. This informs us that most of the hue channel consists of flat, featureless regions. There is a small, but a not-insignificant percentage of pixel gradients  $\in[2,9]$.
This suggests that there are some regions of the hue channel that may represent a hue texture.
Gradients $>9$ appear to be rare in the data set with percentages of pixels $<1\%$.
Due to the decreasing percentage of pixels, it becomes increasingly less likely that these gradients represent a significant region of an image and more likely that they represent the edge of a region of hue. For the multitude of plausible hues, the gradient at these edges will vary significantly over an unknown distribution. For the SHSM measure, we want to treat all edge pixels as equally important. For this reason, we process the gradient values with a logistic function, see Equation (\ref{logistic}). With this, all values below $2$ will be treated as close to zero and values above $9$ will be close to $1$ and therefore all treated the same (see right-hand vertical axis Fig. \ref{textureVBoundary}).  Values in between will be given a sliding scale, which is close to linear, via the $k=1$ and a central value of $h_{0}=5$. These values will be assumed to be a texture that is visible to the viewer and therefore, should be similar in scale between the ground-truth  and a plausible colourisation.

\begin{figure}[!t]
\centering
\includegraphics[width=.5\textwidth]{images/fig1.png}
\caption{Here we see the histogram of occurrences of gradients in hue space for various choices of $c_0$. In the unprocessed case, we see several spikes in the number of occurrences whose positions can all be explained by small values in the a*b* channels. We have highlighted $c_0=5$ as the first of the curves that seem to remove these spikes, but otherwise not further concentrate the number of occurrences in the smaller gradient values.}
\label{removeNuicance}
\end{figure}

\begin{figure}[!t]
\centering
\includegraphics[width=.5\textwidth]{images/fig2.png}
\caption{Here, the black centre line of the curve represents the histogram of gradient values using $h_0=5$, and the grey area represents the standard deviation of the number of occurrences over the 4000 images. We have superimposed the logistic function, and the right-hand y-axis shows what input gradient values will become at the output of the logistic.}
\label{textureVBoundary}

\end{figure}

\begin{figure*}[!t]
            \includegraphics[width=0.22\textwidth]{images/GTP.png}
           \includegraphics[width=0.22\textwidth]{images/CLPA.png}
            \includegraphics[width=0.22\textwidth]{images/CLPB.png}
            \includegraphics[width=0.22\textwidth]{images/CLPC.png}

            \hspace{0.225\textwidth}\includegraphics[width=0.22\textwidth]{images/CLPD.png}
            \includegraphics[width=0.22\textwidth]{images/psGreyL.png}
            \includegraphics[width=0.22\textwidth]{images/L2stripes.png}
    \caption{%
        A ground-truth  image (left) with four plausible colourisations produced in Adobe Photoshop\textsuperscript{TM} and two implausible colourisations.
     }%
    \label{psRecolour}
\end{figure*}

\begin{table*}[!t]
\begin{center}
 \begin{tabular}{|l|l|l|l|l|l|l|}
 \hline
 & \multicolumn{2}{c|}{MSE } &\multicolumn{2}{c|}{PSNR } & \multicolumn{2}{c|}{Proposed Method }\\
 \hline
                 & Hue & Chroma & Hue & Chroma & Hue & Chroma\\
                \hline
Colourisation 1 & 2610     &  394   & 13.96dB    &  22.17dB   & 0.902     & 0.915    \\
\hline
Colourisation 2 & 1660     &  22.41   & 15.9dB      &  34.6dB  & 0.944      &  0.962   \\
\hline
Colourisation 3 & 592     &  110   & 20.41dB     &  27.70dB   & 0.956      &  0.958   \\
\hline
Colourisation 4 &  2485     &   138 & 14.17dB     &  26.74dB   & 0.935      &  0.940   \\
\hline
Stripe Pattern      & 2588     &   728  & 14.00dB      &  19.51dB  & 0.484   & 0.208    \\
\hline
Grey-scale  & 2268     &  2152    & 14.57dB     &   14.8dB  & 0.499      &   0.033 \\
\hline
\end{tabular}
\caption{Four plausible Adobe Photoshop\textsuperscript{TM} colourisations and two implausible colourisations are compared to a ground-truth image. MSE/PSNR fail to separate plausible from implausible colouristions while the proposed method does.}
\label{psTableResults}
\end{center}

\end{table*}

\begin{figure*}[!t]
    \begin{center}
	\includegraphics[width=0.75\textwidth]{images/SSIM_SSHM_map.png}
	\end{center}
	\caption{The left-hand column is the first colourisation and the middle column the second colourisation, there is no ground-truth in this case. We have attempted, with user hints to make a red and green version. The third column shows the similarity using our method, where 0 (black) represents minimum similarity, and 255 (white) represents maximum similarity.  The hue has been penalised, particularly at the edges (middle row). We see that the chroma values produced are not the same (bottom row).}
	\label{mortorPestle}
 \end{figure*}

\begin{table*}[!t]
\begin{center}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
 & \multicolumn{2}{c|}{MSE } &\multicolumn{2}{c|}{PSNR } & \multicolumn{2}{c|}{Proposed Method }\\
 \hline
 Mortor and Pestle image & Hue & Chroma & Hue & Chroma & Hue & Chroma\\
\hline
 Colourisation Comparison & 1634     &  27.21  & 16.0dB    &  33.78dB   & 0.963     &  0.978   \\
\hline

\end{tabular}
\caption{A comparison of two user-guided colourisations using the method of \cite{ZhangZIGLYE17UserGuidedColorization}, see Fig \ref{mortorPestle}}
\label{mortorPestleTableResults}
\end{center}

\end{table*}

\begin{figure*}[!t]

        \includegraphics[width=0.5\textwidth]{images/SOTAhue.png}
        \includegraphics[width=0.5\textwidth]{images/SOTAchroma.png}

	\caption{For 20 images from the Places dataset, we compare six methods' colourisations against the ground-truth .
	ChromaGAN, \cite{Vitoria_2020_WACV}, DeOldify\cite{DeOldify}, Iizuka,  \cite{IizukaLetThereBeColor:2016:LCJ:2897824.2925974}, Larsson, \cite{larsson2016learning}, Nazeri, \cite{nazeri2018image}, Zhang 1, \cite{ZhangIE16ColorfulImageColorization}, Zhang 2,  \cite{ZhangZIGLYE17UserGuidedColorization} (which is used without any user hints in this case).}
	\label{boxPlots}

 \end{figure*}

\section{Experimental Results}
Ideally, we would like to have multiple colour versions of each natural image to test the measurement system. Unfortunately, there is no ideal algorithm for colourisation or re-colourisation. Instead, we will firstly recolour specific sections of a ground-truth image in Adobe Photoshop\textsuperscript{TM} using the magic wand tool to select segments and then change the hue \footnote{Hue in Adobe Photoshop \textsuperscript{TM} will not be from the same space as $L^*h_{ab}^{*}c_{ab}^{*}$, so it will not only affect the hue in our $L^*h_{ab}^{*}c_{ab}^{*}$ experiments.}  of those segments. To the human eye, these look plausible, see Fig. \ref{psRecolour}. We will expect these to score well on similarity. We also include a grey-scale image, and an image with a stripe pattern, with each of the stripes being close to the original ground-truth colour. These look like poor colourisations to the human eye so should score poorly. Table \ref{psTableResults} shows a comparison of the score of each against the ground-truth on MSE, PSNR and our method for hue and chroma channels.
For MSE, a lower number indicates more similarity. For PSNR and our method, the higher number suggests more similarity. While the change made in Adobe Photoshop\textsuperscript{TM} was only to hue, it also changed chroma. We see that both MSE and PSNR have quite a variability in chroma scores for the four colourisations, but they are better than the score for the stripe pattern or grey-scale as would be expected. With hue for MSE and PSNR, we see that some of the Adobe Photoshop\textsuperscript{TM} colourisations receive a worse or on-par score to the stripe and grey-scale.
For our method (standard SSIM) the chroma for the four colourisations has relatively small variance and is on the opposite end of the score range to the grey-scale and stripe pattern which our method has heavily penalised. For hue, we once again see broadly similar scores for the four plausible colourisations. The stripe pattern and grey-scale are much lower, but we note, not close to zero. This is because there are still many parts on the image where the gradient of hue is the same (i.e. small or zero), and this registers as similar. In general, we would advise taking the product of the two values (hue and chroma) as a total score of similarity of two colourisations.

In a second test, we do a qualitative analysis of what our algorithm penalises. We will use the work of  \cite{ZhangZIGLYE17UserGuidedColorization} with their user-guided colourisation, see Fig. \ref{mortorPestle}, to make two colourisations of a grey-scale image. There is a combination of user variability and non-determinism in the system that makes it difficult to produce two colourisations that are close to being as spatially consistent as the system may ultimately be capable of achieving. We want to use our system to see how similar it grades the two colourisations and to show where and why it penalises them. Note: we are comparing these against each other, but the order of the images in the comparison does not matter as the measure is symmetrical.
 Fig. \ref{mortorPestle} shows
 that the chroma (middle row) is not the same for the green and red images and the algorithm has penalised areas in the white stripes in the hue comparison (bottom row), particularly at the borders. A close inspection of the green image shows that the green colour has bled into the white, which has not happened in the red image. Our algorithm has penalised for a section on the right of the top (larger) colour stripe. A close inspection shows a red tint in the green image, which has no correspondence in the red image. Despite this, the two images are very similar, and this does achieve a high score of 0.946 when chroma and hue results are combined by multiplication. Full results for the comparison of the two images in chroma and hue for MSE, PSNR and our method are shown in table \ref{mortorPestleTableResults}. Once again, MSE and PSNR would not rate these images as similar, but our method does.

 Our third test is to see how some of the recent deep-learning methods fare on our measure.
We use a small sample of 20 images\footnote{Some of the methods are only available via a web interface which doesn't allow for automated testing of large sets of images} from the Places dataset to test each method. We compare the models of \cite{DeOldify}, \cite{IizukaLetThereBeColor:2016:LCJ:2897824.2925974}, \cite{larsson2016learning}, \cite{nazeri2018image},\cite{Vitoria_2020_WACV},   \cite{ZhangIE16ColorfulImageColorization}, \cite{ZhangZIGLYE17UserGuidedColorization} (which is used without any user hints in this case).
The results for chroma (SSIM) and hue (SHSM) are shown in Fig. \ref{boxPlots}.

We can see that there is considerable variability in all methods, even over a small sample of 20 images. We should note that the results can be negatively affected by the compression quality of output from the system. DeOldify and Nazeri in particularly showed significant degradation of the L-channel when SSIM was used even though the L-channel should be identical. This is likely caused by compression in JPEG files, where some of the other methods allowed for png file output. The method of  \cite{ZhangIE16ColorfulImageColorization} and \cite{Vitoria_2020_WACV} also used JPEG but the file sizes are larger, suggesting less compression. JPEG compresses the colour information significantly more than the grey-scale information so we can expect that this negatively affects some of the methods' results.

\section{Future Work}
Our simplifying assumption that chroma would share a single-mode between all colourisations, while necessary for the present, is difficult to verify. We hope that through experimentation with this current measure, we can explore and reduce our uncertainty with the chroma channel.

The method we presented here does not attempt to measure whether the hue applied to a semantic segment is plausible for a natural image. This would require knowledge of all object types and the probability distributions of colour tones for each, and measuring the divergence of the predicted distributions from the plausible distributions.
ChromaGAN, \cite{Vitoria_2020_WACV}, attempted to learn this distribution implicitly by minimising the KL divergence of the class distribution with a pre-trained VGG16 network, but this could be attempted explicitly with the hue distribution for classes.

As well as a measure, it may be possible to use our method directly to minimise the loss of a colourisation network to improve spatial coherency in colourisation.

Our work is a full-reference performance measure. Following the assumptions of \cite{LevinColorizationOptimization2004:CUO:1015706.1015780}, \cite{SapiroInpaintingColors1530151} and \cite{YatzivChrominanceBlending1621234}, thought could be given to a spatial comparison between the L-channel prior and the predicted hue, thus enabling a no-reference performance measure.

\section{Conclusions}

We showed that widely-used pixel-difference measures of colourisation plausibility such as MSE and PSNR are misleading.  The hue channel is particularly problematic for these measures, as plausible colourisations can be drawn from a distribution of hues, but converting to other colour spaces propagates the hue problems to other channels. Using an $L^*h_{ab}^{*}c_{ab}^{*}$ formulation of CIEL*a*b*, the problems can be to compartmentalised into the hue channel only. We introduced a measure called SHSM for the hue channel, which we show measures the spatial coherence of hue-channel predictions in colourisation. We presented findings that current state-of-the-art CNN based colourisation methods exhibit high variance in spatial coherence. The SHSM approach will prove useful in progressing the field of colourisation, allowing researchers to assess the spatial coherence of their results both quantitatively and qualitatively.

\bibliographystyle{apalike}


\section{INTRODUCTION}

Our hierarchical motor control system is composed of several areas that are responsible for different functions \cite{merel2019hierarchical}. The high-level areas are in charge of decision making and planning sequence of motion. On the other hand, low-level areas control muscles and regulate both forces and velocities while interacting with dynamic environments. The motor cortex transforms the motion intentions into motor commands, then sends them to the brainstem and the cerebellum. Research studies have shown that patients with cerebellar lesions suffer from clumsy staggering movements \cite{kandel2000principles}.

Hence, various computational models were built around the cerebellar theories on internal model formation for coordination, precision, and accurate timing \cite{porrill2013adaptive}. In particular, the cerebellum seems to act as an inverse model to generate corrective motor commands, as a forward model to enhance sensory predictions, and as a combination of both models \cite{wolpert1998multiple, corchado2019integration}. Considering the effects of cerebellar lesions on both action and perception, it was suggested that the cerebellum is particularly involved in predicting the next sensory states and that the learning process is driven by the error in sensory states \cite{bastian2011moving}. In this paper, we take advantage of the cerebellum capabilities to process the sensory information and predict sensory states to enhance robotic target reaching actions.

Most of the studies consider the cerebellar learning relying on error in joint space \cite{antonietti2019, abadia2019robot, tolu2020cerebellum, Capoleietal2020}. However, in many cases a reference signal in joint space may not exist and thus, learning can be only be carried out based on task-space error.
Cerebellar models were categorized in \cite{van2016neurorobotics}, which they can be either a state-encoder-driven model, functional or cellular model. More attention is given to functional models to approximate functions executed by each cerebellar layer, which makes it more appealing from computer-science perspective \cite{tolu2020cerebellum}. While the cellular models are more challenging, it is not crucial only to verify and develop theories about the learning mechanisms, but to gain insight as well on how neurons act on individual and population basis to achieve such task \cite{abadia2019robot}. Additionally, the spiking nature of neurons makes it possible to run the model on neuromorphic chips, thus allowing less energy consumption and real-time computing \cite{furber2016large, Bairaetal2017}.

\begin{comment}
In \cite{wilson2016cerebellar}, an algorithm was developed based on the adaptive filter model of the cerebellum to control a non-linear actuator. It introduces a functional model that uses sensory error to teach the cerebellar model to compensate for the nonlinearities in the actuator.
\end{comment}
In \cite{Capoleietal2020}, the introduced model is based on the adaptive filter theory. The developed control scheme combines machine learning and some techniques from computational neuroscience for the sensorimotor optimization. However, it does not model the spiking nature of the cerebellum.
In \cite{abadia2019robot}, a compliant controller was built based on a spiking cerebellar model. The cerebellar model relies on joint-space errors based on reference values from a trajectory generator. Thus, it lacks the ability to directly learn from sensory feedback.

In this work, we developed a cellular-level cerebellar controller that guides the motion of a robot based on real-time sensory feedback information. First, a sensorimotor differential map is coarsely built and trained through a motor babbling process. Next, a forward cerebellar internal model is formed to predict and adjust the discrepancies in sensory readings to minimize the deviation from reference targets.

For the first time, a control system that integrates two spiking neural networks (SNNs) to, respectively, emulate the cerebellum's role in anticipatory and adaptive control, and form the differential sensory maps during a vision-guided motion task for a physical robot, is presented. The SNNs incorporate most of the biological neuronal dynamics which include complex and realistic firing patterns based on spikes \cite{maass1997networks}.
Besides, the novelty of the work is that the forward cerebellar predictions are integrated as a Smith predictor \cite{tolu2020cerebellum} and they are supervised by the error in task space.

To show the effectiveness of our approach, we present a detailed experimental study in which  a robot arm performed vision-guided manipulation tasks like reaching a target and following a contour. Results show that the deviation from the desired path is reduced and the execution time is decreased thanks to the cerebellar sensory predictions, which is coherent with experimental results from humans \cite{dewolf2016spiking}.

In the following sections we are describing the developed SNNs and controller (Sec. II), the experimental setup and results (Sec. III), and finally, we are analysing the results and main conclusions (Sec. IV).

\section{METHODS}\label{sec:methods}

A biologically inspired control system shown in Fig.\ref{fig:chart} for voluntary movements is built by integrating two spiking computational models: the differential map ($DM$) and the cerebellar-like controller ($CB$). The sensorimotor map in $DM$ block is formed thanks to the motor babbling process and it receives the modulated sensory input by the $CB$ outcomes.

\subsection{Cerebellar-Based Smith Predictor}

\begin{figure}[b]
\centering
\includegraphics[width=0.6\linewidth]{figures/chart_.pdf}
\caption{The block diagram of the control system. The cerebellum $CB$ acts as a predictor for the next robot state based on $\dot{\theta}_{cmd}$ generated by the differential map $DM$. The sensors $S$ provide the current state of the robot $R$ with a delay $D$, which is compared to the desired state by the Target Generator block $TG$ to generate $\dot{x}_{ref}$.}
\label{fig:chart}
\end{figure}

The cerebellar computational model proposed in this study forms a forward model to predict the next sensory states based on the current robot states and the desired spatial velocity. The forward model allows the designed controller to act as a Smith predictor, which is known to be capable of compensating long time delays (see Fig.\ref{fig:chart}). A dead time is introduced in a control system due to the time needed for sensing, processing of the inputs, computing the control output and actuation, similar to that in a biological system for the sensory signals to travel through the nervous system for processing and generating an adequate motor command which reaches the muscles. In our control system, the dead time ($D$) is introduced by the delay in sensor readings ($S$) and the time needed for the robot ($R$) to take action.

The $DM$ is the differential map \cite{zahra2020differential} needed to generate motor commands in manner similar to the motor cortex, thus acting as an inverse dynamic model, where motor commands ($\dot{\theta}_{cmd}$) are generated based on the current states (the spatial velocity $\dot{x}_s$ and joint angles $\theta$). In robotics, this is usually expressed in the form of an inverse Jacobian matrix $J_{DM}^{\#}(\theta)$ which estimates the required motor command to obtain a target spatial velocity at the end effector $\dot{x}_{ref}$. While considering the delay $\tau$ mentioned earlier it can be expressed as: $\dot{\theta}(t) = J_{DM}^{\#}(\theta)\dot{x}_s(t-\tau)$ given that no correction is applied to make up for the delay introduced. The $CB$ block acts as a forward model to predict the spatial velocity after applying the generated motor commands:

\begin{equation}
    \dot{x}_{cereb}(t) = \tilde{J}(\theta)\dot{\theta}(t)
\end{equation}

where $\dot{x}_s(t-\tau)$ refers to the delayed sensory readings for the actual cartesian velocity of the end effector, and $ \dot{x}_{cereb}(t)$ is the estimated cartesian velocity based on the command $\dot{\theta}(t)$. In case of an optimal controller, the cerebellar-like model would perfectly resemble the forward model of the robot (ie. $\tilde{J}(\theta) = J(\theta)$), where $J(\theta)$ is the Jacobian matrix that describes the real robot; in the studied case, our cerebellar-like model acts as a predictor to correct the differential map estimations. Thus, $CB$ outcomes approximate $J_{DM}$ itself. The difference between the predicted and actual spatial velocity is used for training the cerebellar model, so it minimizes such difference, such that finally $\tilde{J(\theta)} = J_{DM}(\theta)$ when $\dot{x}_{cereb}=\dot{x}_s$ at all the points in the studied space. The estimations can then be interpreted as:
\begin{align}
\begin{split}
     \dot{x}_{cereb}(t-\tau) &= \tilde{J}(\theta)\dot{\theta}     \\
    &= J_{DM}(\theta)(J_{DM}^{\#}(\theta)\dot{x}_s(t-\tau)) \\
    &= \dot{x}_s(t-\tau)
\end{split}
\end{align}

Then, the $CB$ predicts the next sensory state of the robot ($\dot{x}_{cereb}$) which is compared to the corresponding sensory readings ($\dot{x}_{s}$) in the next cycle. The discrepancy in the prediction ($e_{pred}$) is used to train the $CB$ to give better predictions in future cycles:
 \begin{equation}
     e_{pred} = \dot{x}_s - \dot{x}_{cereb}(t-\tau)
 \end{equation}
 Then, $e_{pred}$ is added to the current cycle prediction to make up for the expected error in current readings to obtain $\dot{x}_{pred}$:
 \begin{equation}
     \dot{x}_{pred} = \dot{x}_{cereb}(t) + e_{pred}
 \end{equation}
 The difference between the $\dot{x}_{ref}$ desired spatial velocity and $\dot{x}_{pred}$ shall provide the necessary correction to make up for the error in the approximations of the differential map ($DM$). The $\dot{x}_{ref}$ is updated by the Target Generator block ($TG$) every cycle by comparing the current position $x_s$ to the target position $x_d$:
 \begin{equation}
     \dot{x}_{ref} = \frac{x_d-x_s}{\parallel x_d-x_s\parallel}
 \end{equation}

\begin{comment}
 The estimation is based on the estimation of the Jacobian $\tilde{J}(\theta)$ which is represented by the plasticity in the $PF$ controlling the firing of $PF$ and thus $DCN$ to enhance the predictions of the end effector's velocity.
\end{comment}

\subsection{The Differential Mapping Spiking Neural Network}

\begin{figure}[!t]
\centering
\begin{subfigure}{0.5\columnwidth}
\centering
\includegraphics[width=\columnwidth]{figures/dmsnn_r_l_kj.pdf}
\caption{}
\label{fig:net}
\end{subfigure}%
\begin{subfigure}{0.5\columnwidth}
\centering
\includegraphics[width=1.0\linewidth]{figures/full_network_breif3.pdf}
\caption{}
\label{fig:train_test_breif}
\end{subfigure}
\caption{(a) The layout of the $DM$ block. (b) A schematic diagram for both the training and control phases for the SNN.}
\label{fig:dmsnn}
\end{figure}

A SNN is composed of two layers (an input and output layer) that are connected through all to all plastic synapses to map two correlated spaces and represent the transformation between these spaces. As the map correlates in this study the velocities in joint space -as an output- and in task space -as an input-, thus it acts as a differential map (block $DM$). For a manipulator of $n$ and $m$ degrees of freedom (DOF) in task space and joint space, $n$and $m$ assemblies of neurons are needed, respectively, each encoding one corresponding dimension of the relative space. Hence, the input layer consists of $n+m$ assemblies; $n$ assemblies ($l^{\dot{x}_{1:n}}$) to encode the spatial velocity ($\Dot{x}\in\mathbb R^n$) and $m$ assemblies ($l^{\theta_{1:m}}$) to encode the angular positions of the joints ($\theta\in\mathbb R^m$), while the output layer consists of $m$ assemblies ($l^{\dot{\theta}_{1:m}}$) to encode the angular velocities of the joints ($\Dot{\theta}\in\mathbb R^m$) as shown in \ref{fig:net}. As in this case each dimension of task space is connected to each dimension of joint space, the network resembles a Jacobian-like transformation where $\Dot{x}=J(\theta)\Dot{\theta}$, such that $J_{ij}=\partial x_{i}/\partial \theta_{j}$,
\begin{comment}
\begin{equation}
J(\theta) =
\begin{bmatrix}
\frac{\partial x_{1}}{\partial \theta_{1}} & \cdots & \frac{\partial x_{1}}{\partial \theta_{m}} \\
\vdots & \ddots & \vdots\\
\frac{\partial x_{n}}{\partial \theta_{1}} & \cdots & \frac{\partial x_{n}}{\partial \theta_{m}}
\end{bmatrix}
\end{equation}
\end{comment}
where $J(\theta)\in\mathbb R^{n\times m}$ is the Jacobian matrix.
The input and output neurons are connected through all to all plastic excitatory and inhibitory synapses. Local inter-inhibitory synapses are added within each assembly of the output layer with increasing inhibition from one neuron to distal neurons.
The plastic synapses are modulated during the training phase as indicated in Fig.\ref{fig:train_test_breif}. The internal sensors (proprioception), such as motor encoders, and external sensors (exterioception), such as cameras, provide the data, collected through motor babbling, needed to be introduced to both the input and output layer during the training phase. The encoding of variables is based on the central (preferred) value $\psi_{c}$ assigned to each neuron, with the contribution of the whole assembly utilizing ``population coding'' \cite{amari2003handbook}).

The synapses are modulated accordingly to form the differential map, and the SNN in $DM$ can be used to guide the robot in servoing tasks after several iterations \cite{zahra2020differential}. However, the SNN in $DM$ represents a coarse map which lacks in accuracy and precision and needs to be modulated to be adequate for fine motion control.

\begin{comment}
\begin{table}[!b]
\caption{$DM$ Network Parameters}
\begin{tabular}{|c||*{5}{c|}}
 \hline
 \multicolumn{6}{|c|}{Properties of neurons} \\
 \hline
 \hline
\backslashbox{Layer}{Parameter}&$a$&$b$&$c$&$d$&$\mathcal{N}$\\ \hline
$l^{\theta_{i}}$&0.1&0.2&-65&2&68\\
$l^{\Dot{x_{i}}}$&0.1&0.2&-65&2&68\\
$l^{\Dot{\theta_{i}}}$&0.02&0.15&-55&6&68\\
 \hline
 \hline
 \multicolumn{6}{|c|}{Synaptic connections} \\
 \hline
 \hline
 \backslashbox{Network}{Parameter}& $S$ & $\tau_{1}$ & $\tau_{2}$ & $C_{I}$ & $C_{E}$\\
 \hline
 SNN & 0.05 & 20 & 18 & -4 & 4\\
 \hline
\end{tabular}
\label{table:net_params}
\end{table}
\end{comment}

\subsection{The Proposed Forward Cerebellar Model}

The cerebellar micro-circuit consists of the granular layer, Purkinje layer, and molecular layer. The granular layer contains the granule cells (GC), which make up most of the cells in the brain (around 80\% of the total number of neurons in the brain \cite{herculano2009human}), along with Golgi cells and the mossy fiber ($MF$) axons. The Purkinje layer contains the Purkinje cells ($PC$), which are considered to be the key player with it's distinctive firing pattern for different inputs. The molecular layer contains the parallel fibers ($PF$) intersecting with climbing fibers ($CF$) (carrying signals from inferior olive) and axons of $PC$, along with basket and stellate cells.

In this work, we modeled the cerebellar microcircuit to achieve fine motor control as it provides corrections for the sensory readings introduced to the differential map developed ($DM$), which is somewhat functionally similar to the cerebellum-motor cortex loop \cite{ishikawa2016cerebro}.
Based on the forward model theory \cite{tanaka2019neural}, the $MF$ encode both the sensory inputs and motor commands, which is the output of the $DM$, where each dimension of either the joint space or task space is encoded by an assembly of neurons. The connections between the $MF$ and the granular layer ($GC$) allows for a sparse coding of the robot states. The inferior olive ($IO$) introduces the teaching signal in terms of the task-space error through the $CF$, calculated as the difference between the predicted and the actual spatial velocities, in this case. By using the plasticity of synapses between $GC$ and $PC$, these connections are modulated under supervision of teaching signals from $CF$ to encode the desired correction at the different states of the robot. The $PC$ the inhibits the corresponding Deep Cerebellar Nuclei ($DCN$) to allow for the right value to be decoded from the activity of these neurons.

The $MF$ is divided into two groups, one for encoding values in joint space ($MF_{JS}$) and the other for encoding values in task space ($MF_{TS}$), having total number of assemblies in the whole $MF$ layer equal to $J_{MF}$. The $PC$, $IO$ and $DCN$ consist of $n_{TS}$ groups of neurons and each group has two assemblies of neurons, one for positive and one for negative change of the corresponding degree of freedom (DOF).

\begin{figure}[!b]
\centering
\includegraphics[width=1.0\linewidth]{figures/cereb_mod_3.pdf}
\caption{(a) A schematic diagram for the cerebellum. The $PCs$, $IOs$ and $DCNs$ are divided into two to represent positive and negative changes in the corresponding DOF. }
\label{fig:cereb}
\end{figure}

\begin{table}[!b]
\caption{$CB$ Network Parameters}

\begin{tabular}{ |p{2.2cm}||p{0.75cm}|p{0.75cm}|p{0.75cm}|p{0.75cm}|p{0.75cm}|  }
 \hline
 \multicolumn{6}{|c|}{Properties of neurons} \\
 \hline
 Area & $a$ & $b$ & $c$ & $d$ & $\mathcal{N}$ \\
 \hline
 $MF$ & 0.1 & 0.2 & -65. &2. & 20\\
 $GC$ & 0.02 & 0.25 & -65. &2. & 1000\\

 $PC$ & 1. & 1.5 & -60. &0. & 8\\
 $IO$ & 0.1 & 0.2 & -65. &2. & 8\\
 $DCN$ & 0.05 & 0.1 & -65. &2. & 4\\
 \hline
\end{tabular}

\begin{tabular}{ |p{2.2cm}||p{1.12cm}|p{1.02cm}|p{1.02cm}|p{1.02cm}|  }
 \hline
 \multicolumn{5}{|c|}{Properties of synapses} \\
 \hline
 Projection & Type & P & $\varepsilon_{init}$ & $\varepsilon_{max}$\\
 \hline
 $MF\rightarrow GC$ & Rand & 4 & 1.6 &-\\

 $MF\rightarrow DCN$ & A2A & 1.0 & 1.0 &-\\

 $GC\rightarrow PC$ & Prob & 0.8 & 0.2 &8.0\\
 $IO\rightarrow PC$ & O2O & - & 500.0 &-\\
 $IO\rightarrow DCN$ & O2O & - & 1.7 &-\\
 $PC\rightarrow DCN$ & O2O & - & -5.0 &-\\
 \hline

\end{tabular}
 \footnotetext  'The Type column defines the pattern of connections from layer A to B  in the 'Projection' column with a parameter value 'P'. 'Rand' means that a P random neurons from A connect to one neuron in B. 'Prob' means that there is a probability P for each neuron in A to connect to each neuron in B. 'A2A' means All-to-All connections, and all neurons from A are connected to All neurons in B. 'O2O' means One-to-One connections, and each neuron A is connected to one corresponding neuron in B.
\label{table:cb_params}
\end{table}

The training of the cerebellum starts after the SNN in $DM$ is trained for enough iterations to form a coarse map of the transformation between the task space and motor space.

The input to the $MF$ is encoded based on population coding using Gaussian receptive fields, and the input current to the $i^{th}$ neuron in the $j^{th}$ assembly of MF $\alpha_{MF_{i,j}}(t)$ can be calculated from the following formula:
\begin{equation}
\alpha_{MF_{i,j}}(t) = \exp {\left(\dfrac{-\Vert{\psi_{MF_j} -\psi_{MF_{i,j}}  \Vert^{2}}}{2\sigma_{MF_j}^{2}}\right)}
\label{eq:f_MF}
\end{equation}
where $\psi_{MF_j}$ is the input value to the $j^{th}$ assembly of MF, and $\sigma_{MF_j}$ is calculated based on number of neurons per layer $N_{l}$, and the range of change of the variable to be encoded from $\Psi_{MF_{j_{min}}}$ to $\Psi_{MF_{j_{max}}}$. $\psi_{MF_{i,j}}$ is the central value of the $i^{th}$ neuron in the $j^{th}$ assembly of $MF$.

Each neuron in the $GC$ is connected to one random neuron from each assembly of neurons in the $MF$ layer. The $MF$ connects, as well, through excitatory synapses to the $DCN$ through all to all connections. The $GC$ then projects to the $PC$ through plastic excitatory synapses with a probability $P_{GC-PC}$, and these connections are the $PF$. The parameters are set such that $PC$ synapses are only modulated when neurons in $IO$ are active. The connection between $IO$ and $PC$, is given by $CF$, which are one to one synapses such that neurons of the same group (ie. same DOF) and same direction (ie. positive or negative change) are connected to each other, to ensure only the right $GC-PC$ synapses are modulated. The $IO$ is connected as well through excitatory synapses to $DCN$. The activity of $IO$ neurons can be given by:
\begin{align}
        \alpha_{IO+_j} =
        &\left\{
            \begin{array}{ll}
          \alpha_{IO_{max}} & e_{pred} > 0 \\
          0 & e_{pred}\leq 0
       \end{array}
        \right.
        \\
        \alpha_{IO-_j} =
        &\left\{
            \begin{array}{ll}
          0 & e_{pred}\ge 0 \\
          \alpha_{IO_{max}} & e_{pred} < 0
        \end{array}
        \right.
\end{align}

where $\alpha_{IO+_j}$ and $\alpha_{IO-_j}$ are the average firing rates of the positive and negative IO assemblies of neurons of the $j^{th}$ DOF, respectively, while $\alpha_{IO_{max}}$ is the maximum firing rate of each neuron in the $IO$.

Similarly, the connections between $PC$ and $DCN$ follow the same concept to enable the activation of the right group of neurons. The output from the $DCN$ is then decoded based on the activation of neurons such that:

\begin{equation}
    \dot{x}_{cereb_j} = \frac {\Sigma \alpha_{DCN+_{i,j}} - \Sigma \alpha_{DCN-_{i,j}}}{{\alpha_{DCN_{max}}} * n_{DCN}} *\dot{x}_{max_j}
\end{equation}

where $\dot{x}_{cereb_j}$ is the predicted velocity of the $j^{th}$ DOF, $\alpha_{DCN+_i,j}$ and $\alpha_{DCN-_i,j}$ are the firing rates of the $i_{th}$ neuron in the positive and negative assemblies of the $j^{th}$ DOF, respectively, $n_{DCN}$ the number of DCN neurons in each assembly and $\alpha_{DCN_{max}}$ is the maximum firing rate of each neuron in the $DCN$.

The network parameters are set to allow the $PF$ synapses to change with an adequate pace and reach a certain value after several iterations, such that enough training data is provided. The final synaptic strength of the plastic connections enables $PC$-$DCN$ inhibitions at the proper timing thanks to the teaching signals provided earlier from the $IO$. The neurons in the $PC$ are chosen to act as bistable neurons such that high firing rates can be obtained, yet the firing action can still be controlled by the teaching signal to induce synaptic gating. For this, an inhibiting signal is introduced all the time to the $PC$ as a substitute of the inhibition introduced in the biological model from basket and stellate cells. The strength of the synaptic connections between $IO$ and $PC$ ($CF$) are set to a big value, analogous to having around 300 synapses from each $CF$ to each $PC$ \cite{van2016neurorobotics} to make sure that each spike in $IO$ is enough to evoke spikes in the $PC$, but yet having one synapse instead to save computational power required to represent redundant signals.

\subsection{Neuron Model}

Izhikevich's simple model of spiking neurons \cite{izhikevich2003simple} is adopted for its capability of reproducing various firing patterns while holding a balance between computational cost and biological plausibility \cite{izhikevich2003simple}.
The model is formulated by a set of differential equations:
\begin{align}
\dot{v} & = f(v,u) = 0.04v^2+5v+140-u+I    \label{eq:update_v} \\
\dot{u} & = g(v,u) = a(bv-u)               \label{eq:update_u}
\end{align}
The after-spike membrane potential is reset based on:
\begin{equation}
\label{eq:reset}
\text{if }v\ge30 \text{ mV}, \quad
\text{then } v\leftarrow c,~ u\leftarrow (u+d)
\end{equation}
where $v$ is the membrane potential and $u$ is the membrane recovery variable.
The parameter $a$ determines the decay rate (time constant) for $u$, $b$ controls the sensitivity of $v$ to the sub-threshold fluctuations,
$c$ describes the after-spike membrane potential, and $d$ gives the after-spike value of $u$.
The term $I$ gives the summation of external currents.

\subsection{Synaptic Connections}
The adopted learning rule is the Spike-timing dependent plasticity (STDP), where increase or decrease in the strength of the connections is dependent upon the relative timing of spikes in pre-synaptic and post-synaptic neurons \cite{stdp}.

\begin{comment}
\begin{figure}[t]
\centering
\begin{subfigure}{0.5\columnwidth}
\centering
\includegraphics[width=0.95\linewidth]{figures/stdp_asymm.pdf}
\caption{}
\label{fig:stdp_asym}
\end{subfigure}%
\begin{subfigure}{0.5\columnwidth}
\centering
\includegraphics[width=0.95\linewidth]{figures/stdp_sym.pdf}
\caption{}
\label{fig:stdp_sym}
\end{subfigure}
\caption{(a) Asymmetric and (b) Symmetric learning STDP rules. For asymmetric learning, $S_a$ and $S_b$ are plotted at different values while keeping $\tau_{a}$ and $\tau_{b}$ constant. For the symmetric one, $\tau_{2}$ is plotted at different values while keeping $\tau_{1}$ constant.}
\label{fig:stdp}
\end{figure}
\end{comment}

The anti-symmetric STDP learning rule \cite{stdp_asym_pic} is applied for the plastic synapses in the cerebellum model, and can be expressed as:
\begin{equation}
\Delta {\varepsilon}_{ij} = \left\{
        \begin{array}{ll}
            -S_{a} \exp \left( {-\Delta t}/{\tau_{a}} \right) & \quad \Delta t \leq 0 \\
            \\
             S_{b} \exp \left( {-\Delta t}/{\tau_{b}} \right) & \quad \Delta t > 0
        \end{array}
    \right.
\label{eq:asym_STDP}
\end{equation}
where $S_{a}$ and $S_{b}$ determine the magnitude of the depression and potentiation of synaptic weights, respectively, while $\tau_{a}$ and $\tau_{b}$ decide the time window for occurrence of depression and potentiation, respectively.

For the $DM$, the symmetric model \cite{woodin2003coincident} is adopted, and can be expressed as:
\begin{equation}
\Delta {\varepsilon}_{ij} = S \left( 1-\left( {\Delta t}/{\tau_{1}} \right)^2\right) \exp\left({|\Delta t|}/{\tau_{2} } \right)
\label{eq:sym_STDP}
\end{equation}
where $S$ decides the magnitude of the change in synaptic weights, the ratio between the two variables $\tau_{1}$ and $\tau_{2}$ controls the time window through which potentiation and depression occurs, and $\Delta t$ is the timing difference of spikes at post-synaptic $t_{post}$ and pre-synaptic $t_{pre}$ neurons.

\begin{comment}
In this study, we use $S = 0.05$, $\tau_{1}=20$ms and $\tau_{2}=18$ms.
In this learning model, the change in synapse's weight is controlled by the absolute value of $\Delta t$, but not the sign, i.e. the sequence of firing.
The chosen time window for the pre- and post-synaptic neurons for STDP is 30ms, which means as long as $-30 \leq \Delta t \leq 30$, it will still contribute to the modification of the synaptic strength.
\end{comment}

\begin{figure}[b]
\centering
\begin{subfigure}{0.49\columnwidth}
\includegraphics[width=\linewidth]{figures/servoing_112.pdf}
\caption{}
\label{fig:robot_setup}
\end{subfigure}
\begin{subfigure}{0.49\columnwidth}
\includegraphics[width=\linewidth]{figures/servoing_111.pdf}
\caption{}
\label{fig:servo_compare}
\end{subfigure}
\caption{(a) The experimental setup for the UR3 robot; sensory data is received from the camera and robot encoders. (b) The target reaching with (green) and without (red) the cerebellum in action.}
\label{fig:experiments}
\end{figure}

\section{RESULTS}\label{sec:method_ver}
\subsection{Setup}

To test the developed fully spiking control system, the shoulder and elbow of a UR3 robot are controlled to move in a planar workspace, as shown in Fig.\ref{fig:robot_setup}. A camera provides the current position of the end-effector, which is then compared to the target position to generate the desired spatial velocity vector. On the other hand, the motor encoders provide both the angular positions and velocities necessary to represent the joint space. The end-effector is traced by a colored marker with the fixed camera opposing to the relative workspace. The SNN in $DM$ and $CB$ are built by using the NeMo package \cite{nemo} which runs as well on CUDA enabled GPUs (GeForce GTX 1080Ti in our case).

The motor babbling is run to collect the required data for training the network, such that the limits of both the joint space and workspace are defined. The shoulder motion range is [-110\degree,-30\degree] and the elbow motion range is [60\degree,150\degree]. Training the SNN in $DM$ requires around 3000 iterations which are equivalent to 240s in simulation time and 30s in real time, where an iteration means introducing the inputs for 80ms and the output is decoded based on neural activity during that period. Afterwards, the training of the cerebellum model starts, and the learning duration depends on the task to be achieved. Three tasks are conducted as follows.

\begin{figure}[!b]
\centering
\includegraphics[width=\linewidth]{figures/rand_servo.pdf}
\caption{Random target reaching.}
\label{fig:rand_servo}
\end{figure}

\begin{figure}[!b]
\centering
\begin{subfigure}{\columnwidth}
\centering
\includegraphics[width=0.9\linewidth]{figures/error_pred_.pdf}
\caption{}
\label{fig:pred_error}
\end{subfigure}
\centering
\begin{subfigure}{\columnwidth}
\includegraphics[width=0.9\linewidth]{figures/dcn_firing_.pdf}
\caption{}
\label{fig:dcn_firing}
\end{subfigure}
\caption{(a) Discrepancy in predictions of spatial velocity in both $x$ (${e^{x}}_{pred}$) and $y$ (${e^{y}}_{pred}$) directions. A target is chosen randomly, and each training cycle represents a motion attempt from the starting point to the target. (b) Spike raster of neurons in the DCN layer at the beginning and the end of the training for reaching the target.}
\label{fig:rand_trgt}
\end{figure}

\subsection{Random Target Reaching}\label{sec:target_servo}
The robot is instructed to move from/to random target points in the workspace by only relying on the SNN in the $DM$, ie. before the cerebellum model is activated in the control system. After ca. 10000 training iterations, the same motions are repeated with the active cerebellum model, as shown in Fig. \ref{fig:rand_servo}, with an example of the reaching process in Fig. \ref{fig:servo_compare}. The maximum deviation from the desired straight path is reduced with the active cerebellum applied by a mean of 31\%, and a standard deviation of 18\%, while the reaching time is reduced by ca. 75\%. While repeating the reaching for a randomly chosen target (with each reach denoted as a training cycle), the discrepancy between the sensory prediction and the real position of the end-effector is decreased as shown in \ref{fig:pred_error}. When the $PF$ synapses had been modulated, the activity in the $DCN$ changes as well such that the spikes act collectively to encode the sensory predictions as seen in Fig. \ref{fig:dcn_firing}. While these predictions are provided ahead (one iteration earlier), it allows for the modulation of the reference value sent to the $DM$ before attempting the motion.

\subsection{Radial Reaching}\label{sec:Sequence_servo}

\begin{figure}[!b]
\centering
\begin{subfigure}{0.32\columnwidth}
\centering
\includegraphics[height=0.92\columnwidth]{figures/star_servoing_init_0c.pdf}
\caption{}
\label{fig:star_servo_init}
\end{subfigure}
\begin{subfigure}{0.32\columnwidth}
\centering
\includegraphics[height=0.92\columnwidth]{figures/star_servoing_mid_0c.pdf}
\caption{}
\label{fig:star_servo_mid}
\end{subfigure}
\begin{subfigure}{0.32\columnwidth}
\centering
\includegraphics[height=0.92\columnwidth]{figures/star_servoing_final_0c.pdf}
\caption{}
\label{fig:star_servo_final}
\end{subfigure}
\caption{The robot moves toward eight target reaching points starting from the same point (0,0) with a 45$\degree$ deviation from each target. The performance is shown in (a) before training the cerebellum; (b) after 9600 iterations of training and 4 repetitions of each target reaching; (c) after 19200 iterations of training and 8 repetitions.}
\label{fig:star_servo}
\end{figure}

\begin{figure}[!b]
\centering
\begin{subfigure}{0.49\columnwidth}
\centering
\includegraphics[width=\columnwidth]{figures/infinity_servo1.pdf}
\caption{}
\label{fig:inf_servo_no_cb}
\end{subfigure}
\begin{subfigure}{0.49\columnwidth}
\centering
\includegraphics[width=\columnwidth]{figures/infinity_servo2.pdf}
\caption{}
\label{fig:inf_servo_cb}
\end{subfigure}
\begin{subfigure}{\columnwidth}
\centering
\includegraphics[width=0.8\columnwidth]{figures/inf_servo_err_.pdf}
\caption{}
\label{fig:inf_servo_err}
\end{subfigure}
\caption{Moving the end-effector in an infinity pattern (a) without and (b) with the cerebellum in action. A comparison of deviation (\textit{e}) from the contour and the execution time is plotted in (c).}
\label{fig:infinity_servo}
\end{figure}

The effect of a cerebellar damage was studied in \cite{Fortier2002CEREBELLARAD}, where patients are asked to move at radially arranged targets to check the ability to make smooth and coordinated movements.
Therefore, in this study the robot is instructed to reach target points located in a eight-angled star with an internal angle of 45$\degree$ that are equally 10 cm distant from the central starting point. The robot is required to repeat the reaching motion eight times for each of the 8 targets to give a better performance as shown in Fig. \ref{fig:star_servo}c compared to what obtained after only 4 repetitions (Fig. \ref{fig:star_servo}b) or without the cerebellum in the system (Fig. \ref{fig:star_servo}a).  The maximum deviation is reduced by a mean of 55\% and a standard deviation of
35\%, while the reaching time is reduced by ca. 120\%. It shall be noted that each target is tested independently. The average number of iterations needed to complete a target reach is around 300, which means that approximately 2400 iterations only are needed to train the network to reach in a specific direction compared to 10,000 iterations of random reaching in the previous subsection.

\subsection{Contour Following}\label{sec:track_servo}
The robot end-effector follows a predefined reference contour, an infinitive- or eight-shape ($\infty$) contour (see Fig.\ref{fig:infinity_servo}) defined by $x_{t} = 0.5R_{t}\sin(2\gamma)$ and $y_{t} = R_{t}\cos(\gamma)$, where $R_{t}$ is the radius of each half of the $\infty$-shape (set to 7cm in this experiment), and $\gamma \in [0,2\pi]$ divided over equal intervals to dicretize the contour into a set of points (80 points in this case) as sequential targets. The end-effector shall attempt to reach the point with a defined accuracy (less than 1mm) within a time limit (1$s$), otherwise the target is skipped and the end-effector moves to the next one. The results of the same task carried out by only the $DM$ and with also the $CB$ in action can be seen in Fig.\ref{fig:inf_servo_no_cb} and Fig.\ref{fig:inf_servo_cb}, respectively. The error is defined in this case as the distance between the end-effector position and the closest point on the contour while moving, and it is shown in the plot in Fig.\ref{fig:inf_servo_err} after filtering by a first order filter. It is clear from the plot that the reduction in both the error and time needed to complete the contour following task starting from one point and ending back at the same point. The maximum error is reduced from approximately 8mm to 2.5mm, and the time to complete the contour is reduced from approximately 46$s$ to 33$s$.

\section{CONCLUSIONS}\label{sec:conclusions}
In this study, we suggested a biomimetic control system for vision-based control that integrates a cellular-level spiking cerebellar model that helps to obtain more accurate and faster robot motions. The teaching signal for the cerebellum is based on sensory feedback (in task-space), and thus it allows task-based adaptation in contrast with the previous literature. The sensory corrections are then fed into a motor-cortex-like spiking network which acts as a differential map to convert the motion plan (in spatial coordinates) into motor commands in joint-space. In the developed control system, the cerebellum acts as a Smith predictor to make up for the lack in accuracy in the developed differential map.

From the results, we conclude that both the accuracy (shown as deviation from the desired path) and execution time are reduced thanks to the sensory predictions of the cerebellum. The performance of the robot in reaching a specific target and reaching random targets indicates that the cerebellar-model is able to learn fast (within only 8 repetitions) a specific example and to also generalize its learning to improve the performance across the whole workspace from random motions. Moreover, from the performed experiments, it shall be noted that some regions in the $DM$ have higher error in reaching than other regions. This is mostly related to the randomness in the followed motor babbling technique. Hence, even while providing sensory corrections from the $CB$, the $DM$ still does not hold much improvement at the correct sensory values. This urges the need to apply/develop a different motor babbling algorithm.
The experimental results obtained with the robot arm display a fair correlation with biological observations \cite{Fortier2002CEREBELLARAD, dewolf2016spiking}. Future aim will be to test the effect of cerebellar disfunctions in the control tasks.
Additionally, the current study only considers planar motion of the robot in line with experiments carried out on patients with cerebellar damage. In the future, non-planar motions needs to be studied along with an algorithm for self-tuning the parameters of the cerebellar model.
While the developed cerebellar model copies just some of the features of the real cerebellum, more biological features are to be included in the future to exploit the full potential of cellular-level models.