\section{Introduction}
\label{sec:intro}
Human pose estimation is an important computer vision problem \cite{johnson2010clustered}. Its basic task is to find the posture of a person via recognising human joints and rigid parts from normal RGB images. The extracted pose information is essential to modelling and understanding the human behaviours, and can be used in many vision application problems, such as virtual/augmented reality, human-computer interaction, action recognition and smart perception.

In the psst few years, pose estimation methods based on deep neural network techniques have achieved great progress\cite{wei2016convolutional}\cite{newell2016stacked}\cite{luvizon20182d}. Although the performance of some human pose estimation models (e.g. \cite{chu2016structured}\cite{chu2017multi}\cite{sun2019deep}) is almost saturated on the above mentioned datasets, applying these high-precision algorithms to the other specific industrial tasks shows a degradation in accuracy. For instance, one application case is workouts or sports scoring. In this case, lots of activities have severe self-occlusion or unusual postures, such as push-up and crunch. We find out the models \cite{he2017mask}\cite{cao2017realtime}\cite{osokin2018real} trained on the MS-COCO dataset \cite{lin2014microsoft} cannot correctly detect body joints with atypical postures, as shown in Fig.~\ref{falseSamples}. In the top-right image of Fig.~\ref{falseSamples}, the right knee is falsely detected as left knee. In the top-left and lower-part images of Fig.~\ref{falseSamples}, some body joints, such as shoulders, knees and ankles, are missed in prediction. Since the pose estimation results of the same person in the standing posture are correct, we argue the false predictions are caused by the abnormal postures. Current datasets lack the corresponding samples\cite{andriluka20142d}\cite{andriluka2018posetrack}.

\begin{figure}[t]
\centering
\includegraphics[scale=1]{false_samples.pdf}
\caption{Limitations of applying current pose estimation models on some workout postures, which have severe self-occlusion. Some body keypoints are falsely detected or missed in prediction even the background is plain.}
\label{falseSamples}
\end{figure}

We use the MS-COCO dataset \cite{lin2014microsoft} as an example to analyse the distribution of human postures. In our statistics, the number of human instances in standing posture achieves 102,495 (84.53\%) while people in other postures only have 18,756 (15.47\%) as shown in Fig. \ref{figSampleDistribution}. The human instances in a horizontal position or an uncommon pose are extremely rare. This makes the model unable to learn the knowledge of irregular postures during training.

\begin{figure}[t]
\centering
\includegraphics[scale=0.2]{SampleDistribution.pdf}
\caption{The posture distribution of MS-COCO dataset. Around 85\% human instances are standing with good, upright posture.}
\label{figSampleDistribution}
\end{figure}

To improve the performance of human pose estimation in the certain sports situation, a large-scale human keypoints benchmark is presented in this paper. Our benchmark significantly advances state-of-the-art in terms of particular activities, and includes more than 16,000 images of people. We used the workout class videos as a data source and collected images and image sequences using queries based on the descriptions of more than 200 workout activity types. For each activity type, there are 3 different shot angles. This results in a diverse set of images covering not only different workout activities, but contrasting postures. This allows us to enhance the current human pose estimation methods.

\section{Related Work}
There are several human keypoints datasets presented in the past decades. Buffy dataset \cite{ferrari2008progressive} and PASCAL stickmen dataset \cite{eichner2009better} only contain upper-bodies, but we need to process the full-body. In these two datasets pose variation is insignificant. The contrast of image frames is relatively low in the Buffy dataset.

The UIUC people dataset \cite{tran2010improved} contains 593 images (346 for training, 247 for testing). Most people in the images are playing badminton. Some people are playing jogging, Frisbee, standing, walking, etc. There are very aggressive pose and spatial variations. However, the activity type is limited in this dataset.

The sport categories of Sport image dataset \cite{wang2011learning} is more plentiful, which including soccer, cycling, acrobatics, American football, croquet, golf, horseback riding, hockey, figure skating, etc. The total number of images is 1299 (649 of them are split as training set and the rest as testing set).

Leeds Sports Poses (LSP) dataset \cite{johnson2010clustered} includes 2000 images, where one half for training and the other half for testing. The dataset shows people involved in various sports.

The image parsing (IP) dataset \cite{ramanan2006learning} is a small dataset and contains 305 images of fully visible people, where 100 images for training and 205 images for testing. The dataset consists of various activities such as dancing, sports and acrobatics.

The MPII Human pose dataset \cite{andriluka20142d} consists of 24,589 images, in which 17,408 images with 28,883 annotated people are split for training. During the testing stage, one image may contain multiple different evaluation regions that consist of a non-identical number of people. \cite{pishchulin2016deepcut} defines a set of 1,758 evaluation regions on the test images with rough position and scale information. The evaluation metric deploys mean Average Precision (mAP) of the whole body joint prediction. The accuracy results are evaluated and returned by the staff members of the MPII dataset.

The MS-COCO keypoints dataset \cite{lin2014microsoft} includes training, validation and testing sets. On the COCO 2017 keypoints challenge, training and validation sets have 118,287 and 5000 images respectively, totally containing over 150,000 people with around 1.7 million labelled keypoints. In experiments, we perform ablation studies on the validation set. To analyse the effect of training, we also combine the COCO train set with the FollowMeUp train set to validate that new images will not affect the model's generality performance.

The DensePose-COCO dataset \cite{alp2018densepose} has reannotated dense body surface annotations on the 50k COCO images. These dense body surface annotations can be understood as continuous part labels of each human body.

The PoseTrack dataset \cite{andriluka2018posetrack} includes both multi-person pose estimation and tracking annotations in videos. It can perform not only pose estimation in single frames, but also temporal tracking across frames. The dataset contains 514 videos including 66,374 frames in total. The annotation format defined 15 body keypoints. For the single-frame pose estimation, the evaluation metric uses mean average precision (mAP) as is done in \cite{pishchulin2016deepcut}.

\section{The Dataset}

\subsection{Pose Estimation}
The key motivation directing our data selection strategy is that we want to represent rare human postures that might be not easily accessed or captured. To this end, we follow the method of \cite{ainsworth20112011} to propose a two-level hierarchy of workout activities to guide the collection process. This hierarchy was designed according to the body part to be trained during the exercise. The first level is the body part interested to be trained, such as shoulder, whereas the second level is specific workout activities that can strengthen the muscles of shoulder.

\subsubsection{Data collection} We select candidate workout videos according to the hierarchy and filter out videos of low quality and those that people are truncated. This resulted in over 600 videos spanning over 200 different workout types with three shot angles. We also filter out the frames in which pose is not recognisable due to poor image quality, small scale and dense crowds. This step resulted to a total of 110,000 extracted frames from all collected videos. Secondly, since different exercises have disparate periods, we manually pick key frames with people from each video. We aim to select frames that either depict the whole one exercise period in a substantially different pose or different people with dissimilar appearance. The repeated or no significant distinction postures are ignored. Following this step we annotate 16,519 images. We rough randomly split the annotated images for training and use the rest for testing. Images from the same video are either all in the training or all in the test set. We last obtain the train set of 15,435 images and test set of 1,084 images.

\subsubsection{Data annotation} We follow the keypoint annotation format of COCO dataset, where 17 body keypoints are defined. This design facilitates us to utilise the common samples of COCO dataset during training. Following \cite{lin2014microsoft} the left/right joints in the annotations refer to the left/right limbs of the person. Additionally, for all body joints the corresponding visibility is annotated. At test time both the accuracy of joints localisation of a person along with the correct match to the left/right limbs are evaluated. The annotations are performed by in-house workers and inspected by authors. For some unqualified and incorrect annotations are modified continuously until totally correct. To maintain the quality of annotations, we arranged a number of annotation training classes for all annotation workers to unify the standard of annotation. We also supervise and handle some uncertain cases for workers during annotation.

\subsubsection{Pose Estimation Evaluation Metrics} Some previous keypoints evaluation metrics rely on the calculation of body limbs' length, such as ${\rm PCP}$, ${\rm PCK}$ and ${\rm PCK_h}$ used in \cite{andriluka20142d}. However, the workout activities usually have specific postures where the limb's length may be near 0 if the limb is perpendicular to the image plane and the evaluation is not numeric stable in these cases. Therefore comparing the distance between points of groundtruth and prediction directly is more sensible. Here we follow the COCO keypoints dataset, using 5 metrics to describe the performance of a model. They are AP (i.e. average precision), ${\rm AP^0{}^.{}^5}$, ${\rm AP^0{}^.{}^7{}^5}$, ${\rm AP^M}$, ${\rm AP^L}$, as illustrated in Table \ref{tab:cocoMetric}. In the matching between predictions to groundtruth, a matching criterion called object keypoint similarity (OKS) is defined to compute the overlapping ratio between groundtruth and predictions in terms of point distribution \cite{lin2014microsoft}. If OKS is larger than one threshold value (e.g. 0.5), the corresponding groundtruth and prediction are considered as a matching pair and the correctness of predicted keypoint types is further analysed. Here OKS is similar to the intersection over union (IoU) in the case of object detection. Thresholding the OKS adjusts the matching criterion. Notice that in general applications, ${\rm AP^0{}^.{}^5}$ gives a good accuracy already. When computing AP (averaged across all 10 OKS thresholds), 6 thresholds exceed 0.70 are over strict due to unavoidable jittering in annotations.

\begin{table}[!t]
\caption{\label{tab:cocoMetric}Evaluation metrics on the COCO dataset.}
\centering
\begin{tabular}{|p{1.2cm}|p{7.2cm}|}
\hline
Metric & Description\\
\hline
AP & AP at OKS$^*$ = 0.50 : 0.05 : 0.95 (primary metric)\\
\hline
${\rm AP^0{}^.{}^5}$ & AP at OKS = 0.50\\
\hline
${\rm AP^0{}^.{}^7{}^5}$ & AP at OKS = 0.75\\
\hline
${\rm AP^M}$ & AP for medium objects: $32^2 < area < 96^2$\\
\hline
${\rm AP^L}$ & AP for large objects: $area > 96^2$\\
\hline
\multicolumn{2}{@{}l}{$^*$OKS--Object Keypoint Similarity, same role as IoU}
\end{tabular}
\end{table}

\section{Analysis of The State of The Art}
In this section we first compare the leading human pose estimation methods on the COCO keypoints dataset, and then analyse the performance of these approaches on our benchmark.

The basis of the comparison is that we note that there is no uniform evaluation protocol to measure the performance of existing methods from a view of practical application. Although human pose estimation is one of the longest-lasting topics, and significant performance improvement has been achieved in the past few years, some reported accuracies in these approaches are obtained through several post-processing steps or some strategies used in the dataset challenge. For example, performing multi-scale evaluation, refining results by a different method, or precision is evaluated at one image scale while speed is recorded at another scale. These post-processing steps interfere the judgement in identifying the strength and weakness of an algorithm. Therefore, evaluating a method without any post-processing steps and strategies is more objective and more valuable for the research and practical application.

The aim of the analysis is to evaluate the generality of the current models on the different datasets and their performance to the unseen samples, identify the existing limitations and stimulate further research advances.

Currently, there are two main categories of solutions: top-down methods \cite{sun2019deep}\cite{fang2017rmpe}\cite{papandreou2017towards}\cite{xiao2018simple}\cite{chen2018cascaded}\cite{su2019multi} and bottom-up methods \cite{cao2017realtime}\cite{osokin2018real}\cite{insafutdinov2016deepercut}\cite{newell2017associative}\cite{papandreou2018personlab}\cite{kreiss2019pifpaf}. Top-down methods can be seen as a two-stage pipeline from global (i.e. the bounding box) to local (i.e. joints). The first stage is to perform human detection and to obtain their respective bounding boxes in the image. The second stage is to perform single person pose estimation for each of the obtained human regions. \cite{sun2019deep} deploys multiple high-to-low resolution subnetworks with repeated information exchange across multi-resolution subnetworks. This design obtains rich high-resolution representations and leading more accurate result. \cite{fang2017rmpe} utilises a Symmetric Spatial Transformer Network to handle inaccurate bounding boxes. \cite{xiao2018simple} uses simple deconvolution layers to obtain high-resolution heatmaps for human pose estimation. On the side of bottom-up methods, \cite{cao2017realtime} proposes a limb descriptor and an efficient bottom-up grouping approach to associate neighbouring joints. \cite{osokin2018real} modifies the network architecture of \cite{cao2017realtime} and optimises the post-processing steps to achieve real-time speed on the CPU devices. \cite{kreiss2019pifpaf} designs two new descriptors based on \cite{cao2017realtime} for body joints and limbs with the additional variable of object's spread. \cite{newell2017associative} presents a network to simultaneously output keypoint detections and the corresponding keypoint group assignments. \cite{moon2018posefix} designs a feedback architecture that combining the keypoint results of other pose estimation methods with the original image as the new input to the human pose estimation network. In our analysis we consider 8 state-of-the-art multi-person pose estimation methods, which are listed in Table \ref{tab:compare_coco}.

We compare the performance of each approach in terms of accuracy and speed on the COCO dataset and our novel FollowMeUp dataset. All the experiments are performed on a desktop with one NVIDIA GeForce GTX-2080Ti GPU. Since all testing approaches are trained and optimised on the COCO dataset, their open source codes have the corresponding configurations, we directly use their default parameters in our testing.

\subsection{Comparisons of Approaches on the COCO Dataset}
Table \ref{tab:compare_coco} presents the comparison results of testing approaches on the COCO dataset. The upper part of Table \ref{tab:compare_coco} are top-down approaches. \cite{sun2019deep} has the highest AP precision of 0.753. Note that the runtime costs around 50 ms as this only includes the part of pose estimation since this open source library uses the groundtruth of human bounding box as the human detection results on the COCO validation set. \cite{xiao2018simple} and \cite{fang2017rmpe} have a relatively lower accuracy than \cite{sun2019deep} using smaller input sizes, which illustrates that the high-resolution and detailed representation is important for the task of human pose estimation. Note that some post-processing strategies, such as multi-scale and flip, are ignored to obtain the actual performance in the real application environments.

\begin{figure}[t]
\centering
\includegraphics[scale=0.6]{predAnalysis.pdf}
\caption{The comparison of the numbers of effective instance predictions and body keypoints between top-down and bottom-up methods. The prediction number of top-down method is around 10 times higher than bottom-up method.}
\label{figPredAnalysis}
\end{figure}

For the bottom-up methods, \cite{cao2017realtime} achieves the fastest speed. \cite{kreiss2019pifpaf} attains the highest precision in this group. The joint grouping part of \cite{kreiss2019pifpaf} costs much longer time than \cite{cao2017realtime}. \cite{osokin2018real} has around 7\% degradation compared with \cite{cao2017realtime} due to using a light-weight network architecture. We also see that the precision of bottom-up algorithms are lower than top-down methods. After detailed analysis, we find that the numbers of predicted effective keypoints of bottom-up methods are around 10 times less than top-down methods as illustrated in Fig. \ref{figPredAnalysis}.
We note that top-down methods correspond to performing single-person pose estimation on each detected human region. Single-person pose estimation can output all types of keypoints even the keypoint is occluded or truncated. However, for multi-person bottom-up methods, two or more overlapping keypoints with the same type can only be detected one due to depth information is not available on the RGB image. For the COCO dataset, there are a lot of crowded and occluded human instances. Therefore, the performance of bottom-up methods is weakened. In the FollowMeUp dataset, the crowding case is rare while most human instances have self-occlusion. We perform the same comparison on the FollowMeUp dataset and validate that bottom-up methods have comparable performance to top-down approaches in this circumstance.

\begin{table*}[!t]
\caption{\label{tab:compare_coco}Comparisons of pose estimation results on the COCO 2017 validation set.}
\centering
\begin{tabular}{|p{1.6cm}|p{1.3cm}|p{0.8cm}|p{0.9cm}|p{1cm}|p{0.8cm}|p{0.8cm}|p{1.5cm}|p{1.3cm}|}
\hline
Type & Method & AP & ${\rm AP^0{}^.{}^5}$ & ${\rm AP^0{}^.{}^7{}^5}$ & ${\rm AP^M}$ & ${\rm AP^L}$ & Input Size & Runtime \\
\hline
Top-down & HRNet\cite{sun2019deep} & \textbf{0.753} & \textbf{0.925} & \textbf{0.825} & \textbf{0.723} & \textbf{0.803} & 384x288 & 0.049$^*$ \\
 & Xiao\cite{xiao2018simple} & 0.723 & 0.915 & 0.803 & 0.695 & 0.768 & 256x192 & 0.110 \\
 & RMPE\cite{fang2017rmpe} & 0.735 & 0.887 & 0.802 & 0.693 & 0.799 & 320x256 & 0.298 \\
\hline
Bottom-up & PAF\cite{cao2017realtime} & 0.469 & 0.737 & 0.493 & 0.403 & 0.561 & 432x368 & \textbf{0.081} \\
 & Osokin\cite{osokin2018real} & 0.400 & 0.659 & 0.407 & 0.338 & 0.494 & 368x368 & 0.481\\
 & PifPaf\cite{kreiss2019pifpaf} & 0.630 & 0.855 & 0.691 & 0.603 & 0.677 & 401x401 & 0.202\\
 & AE\cite{newell2017associative} & 0.566 & 0.818 & 0.618 & 0.498 & 0.670 & \textbf{512x512} & 0.260\\
 & PoseFix\cite{moon2018posefix} & 0.411 & 0.647 & 0.412 & 0.303 & 0.559 & 384x288 & 0.250 \\
\hline
\multicolumn{2}{@{}l}{$^*$: without human detection}
\end{tabular}
\end{table*}

\subsection{Comparisons of Approaches on the FollowMeUp Dataset}
Table \ref{tab:compare_followmeup} provides the comparison results of testing approaches on the COCO dataset. Since the open source libraries of \cite{sun2019deep} and \cite{xiao2018simple} do not provide default human detection algorithm, using different human detector may bias the precision distribution, thus we do not test \cite{sun2019deep} and \cite{xiao2018simple} on the FollowMeUp dataset. We are surprised that \cite{fang2017rmpe} obtains a very high precision value. However, the training set only including the COCO dataset of \cite{cao2017realtime} just achieve the precision of 0.778. We argue that the training set of \cite{fang2017rmpe} may include other samples except the COCO dataset with particular postures. In this dataset, the precision of \cite{osokin2018real} decreases by 13\% in ${\rm AP^0{}^.{}^5}$ compared with \cite{cao2017realtime}, which indicates that the generality of \cite{osokin2018real} is also narrowed. We use the results of \cite{cao2017realtime} as the initial poses of \cite{moon2018posefix}. Through pose refinement, \cite{moon2018posefix} improved the pose estimation results by 0.4\%.

\begin{table*}[!t]
\caption{\label{tab:compare_followmeup}Comparisons of pose estimation results on the FollowMeUp dataset.}
\centering
\begin{tabular}{|p{1.6cm}|p{1.7cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|}
\hline
Type & Method & ${\rm AP^0{}^.{}^5}$ & ${\rm AP^0{}^.{}^6}$ & ${\rm AP^0{}^.{}^7}$ & ${\rm AP^0{}^.{}^8}$ & ${\rm AP^0{}^.{}^9}$ \\
\hline
Top-down & RMPE\cite{fang2017rmpe} & 0.975 & 0.948 & 0.885 & 0.787 & 0.421 \\
\hline
Bottom-up & PAF\cite{cao2017realtime} & 0.778 & 0.728 & 0.625 & 0.474 & 0.326 \\
 & Osokin\cite{osokin2018real} & 0.645 & 0.585 & 0.520 & 0.370 & 0.215 \\
 & PoseFix\cite{moon2018posefix} & 0.782 & 0.716 & 0.621 & 0.466 & 0.334 \\
\hline
\end{tabular}
\end{table*}

\subsection{The Effect of Training on the FollowMeUp Dataset}
To validate the effectiveness of samples with particular postures, we retrain the model on the COCO + FollowMeUp train set using the method of \cite{cao2017realtime}. Testing is performed both on the FollowMeUp test set and COCO validation set. The results of testing are provided in Table \ref{tab:train_followmeup}. We notice that the performance of the retrained model is greatly improved by around 20\% in ${\rm AP^0{}^.{}^5}$. While the threshold of AP becomes more strict, the AP value is decreased. Even in the most strict threshold of 0.9, the AP value attains 0.691, which is higher than the model before retraining by 37\%. The accuracy comparison of before and after retraining on the FollowMeUp dataset is shown in Fig.~\ref{ap_compare}. We also perform testing on the COCO validation set using before and after retraining models to check whether the model can maintain the performance on the COCO dataset. In Table \ref{tab:test_followmeup} we see that before and after retraining the precision has no change. The generality of the retrained model is preserved. These results show that increasing some unusual samples which had not been learnt by the model before is an effective way to improve the accuracy in some specific scenes.

\begin{table*}[!t]
\caption{\label{tab:train_followmeup}Comparisons of pose estimation results on the FollowMeUp dataset.}
\centering
\begin{tabular}{|p{1.3cm}|p{3.2cm}|p{1.8cm}|p{0.8cm}|p{0.8cm}|p{0.9cm}|p{0.8cm}|p{0.8cm}|}
\hline
Method & Train Set & Test Set & ${\rm AP^0{}^.{}^5}$ & ${\rm AP^0{}^.{}^6}$ & ${\rm AP^0{}^.{}^7}$ & ${\rm AP^0{}^.{}^8}$ & ${\rm AP^0{}^.{}^9}$ \\
\hline
PAF\cite{cao2017realtime} & COCO & FollowMeUp & 0.778 & 0.728 & 0.625 & 0.474 & 0.326 \\
PAF\cite{cao2017realtime} & COCO + FollowMeUp & FollowMeUp & 0.964 & 0.959 & 0.926 & 0.876 & 0.691 \\
\hline
\end{tabular}
\end{table*}

\begin{figure}[t]
\centering
\includegraphics[width=9.5cm]{ap_compare.pdf}
\caption{Comparison of estimation accuracy before and after retraining on the FollowMeUp dataset. The accuracy of retrained model (marked as green triangles) has an obvious improvement.}
\label{ap_compare}
\end{figure}

\begin{table*}[!t]
\caption{\label{tab:test_followmeup}Comparisons of pose estimation results on the COCO dataset.}
\centering
\begin{tabular}{|p{1.3cm}|p{3.2cm}|p{1.2cm}|p{0.8cm}|p{0.8cm}|p{1cm}|p{0.8cm}|p{0.8cm}|}
\hline
Method & Train Set & Test Set & AP & ${\rm AP^0{}^.{}^5}$ & ${\rm AP^0{}^.{}^7{}^5}$ & ${\rm AP^M}$ & ${\rm AP^L}$ \\
\hline
PAF\cite{cao2017realtime} & COCO & COCO & 0.465 & 0.740 & 0.447 & 0.379 & 0.597 \\
PAF\cite{cao2017realtime} & COCO + FollowMeUp & COCO & 0.465 & 0.748 & 0.454 & 0.373 & 0.605 \\
\hline
\end{tabular}
\end{table*}

\section{Conclusion}
The problem of human pose estimation has obtained a great progress in recent years. This progress cannot be done without the development of large-scale human pose datasets. However, the existing human pose datasets are not sufficient for some particular application environments. In this paper, we propose a new large-scale workout activity human pose dataset, which provides a wide variety of sport exercise postures. We select 8 state-of-the-art multi-person pose estimation approaches and compare their performance on both the popular COCO keypoints dataset and our FollowMeUp dataset. The comparison results show that most methods trained on the COCO dataset do not have ideal performance on the FollowMeUp dataset. We also test the generality of the model using the data of both COCO and FollowMeUp datasets. The test results show that training on the data of both COCO and FollowMeUp datasets will not affect the performance of the model on the COCO dataset but the performance of the model on the FollowMeUp dataset is greatly improved. In the future, we will continue investigate pose tracking\cite{raaj2019efficient}, multi-view action recognition\cite{zhao2019semantic}, and light-weight network design\cite{zf2019fastHumanPose} approaches on the FollowMeUp dataset.


\section{Introduction}

A cellular automaton (CA) model consists of an $n$-th dimensional space partitioned into a discrete subset of $n$-dimensional volumes, which are called cells and are defined in a discrete time. A finite list of possible states is defined for each cell, and each cell has one state. A local neighbourhood is defined for each cell at every time step. The state of a cell can be changed by a \emph{state change rule}, which is a rule that allows the computation of the new state for the cell, and is dependent on other cells in the local neighbourhood \cite{Janssens2010}. Typically, this rule is fixed, that is, it is the same rule for all cells. It does not change over time and it is applied to all cells simultaneously. However, this rule can be stochastic, which means that the new states are chosen according to some probability distribution. This rule can also be applied to each individual cell independently and so the new state of a cell affects the calculation of states in neighbouring cells. Chopard \emph{et al.}$\;$ \cite{Chopard1989} presented the first application of cellular automaton to model diffusion on lattices. Subsequently many other CA approaches were applied to reaction-diffusion problems, see Boon \emph{et al.}$\;$ \cite{Boon1996} and Weimar \cite{Weimar1997}. Then, it seems, there is a pause in the literature on using cellular automata approaches to work with the diffusion equation. Moreover, this type of modelling does not seem to have been applied to work with coupled, nonlinear diffusion problems, which are the focus of this paper.

The \emph{Kirkendall effect} is the name given to the physical phenomenon whereby atomic diffusion occurs via a vacancy exchange mechanism instead of a substitutional or ring mechanism. Recently the Kirkendall effect has been used to create hollow nanostructures, which can be used in a variety of applications. In Aldinger \cite{Aldinger1974} they present the first use of the Kirkendall effect to create hollow structures. The hollow permits their use in transporting drugs and biomolecules and then releasing them in a controlled manner \cite{An2009}. Na \emph{et al.}$\;$ \cite{Na2009} present a review on nanostructrues and magnetic resonance imaging (MRI). These structures have also been proposed to enhance the rate capability and cycling stability in lithium-ion batteries \cite{Wang2012}. Hollow nanoparticles have also been reported to be good catalysts \cite{Kim2002,Li2006}.

In an attempt to understand, and so better control the growth of hollow nanoparticles Ribera \emph{et al.}$\;$ \cite{Ribera2017} rigorously derived governing equations for the substitutional binary diffusion problem. Moreover, under sensible assumptions they reduce these governing equations in order to provide an analytically tractable problem. As a starting point they examine the one-dimensional problem of an insulated bar. In this paper we investigate the same problem but from the cellular automaton standpoint. This will help understand further the physical mechanisms behind the Kirkendall effect.

In the following section we present the model of Ribera \emph{et al.}$\;$ \cite{Ribera2017}, which will be the starting point of the cellular automaton model, presented in Section \ref{sec-CA}. In the limit where the number of cells is large we prove that the CA model reduces to a special case of binary diffusion, where one species diffuses much faster than the other. In the Results section we verify this by showing that at large time and with a fast diffuser the CA model coincides with the continuum model. Further, since the computer speed reduces with $N$ it is clear that the CA model is particularly useful at small scales, such as with nanostructures.

\section{Continuum model for substitutional binary diffusion}
\label{sec-model-reduced}

Let us consider a binary crystalline solid composed of three species: atomic species A, atomic species B, and vacancies V. We label the fast diffuser as species A, and the slow one B. Since we are considering a perfect lattice, that means that the sum of all the fluxes is zero and then it is only necessary to work with the evolution of two species to fully define the problem.

Now consider an insulated one-dimensional bar of length $2l$. At $t=0$ the side $x \in [-l,0]$ is made of material A (and a proportion of vacancies), and the side $x \in [0,l]$ is made of material B (and a proportion of vacancies).

For $t >0$ the diffusion of species is defined by
\begin{align}
\label{gov-eq-1-mathmodel}
\frac{\partial X_A}{\partial t} &= \frac{\partial}{\partial x}\left(D_{AA}^V \frac{\partial X_A}{\partial x} \right) - \frac{\partial}{\partial x}\left(D_{AV} \frac{\partial X_V}{\partial x} \right),\\
\label{gov-eq-2-mathmodel}
\frac{\partial X_V}{\partial t} &= - \frac{\partial}{\partial x}\left(D_{VA} \frac{\partial X_A}{\partial x} \right) + \frac{\partial}{\partial x}\left(D_{VV} \frac{\partial X_V}{\partial x} \right),
\end{align}
where $X_i$ are the mole fractions corresponding to the $i$-th species and the diffusion coefficients $D_{AA}^V$, $D_{AV}$, $D_{VA}$ and $D_{VV}$ vary nonlinearly with $X_i$ \cite{Ribera2017}. In the limit where A diffuses much faster than B,
\begin{equation}
\label{reduced-diffusion}
\begin{split}
\;\;\;\;\;\;\;\: \hat{D}_{AA}^V &\sim \Gamma X_V, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\:\! \hat{D}_{AV} \sim \Gamma X_A,\\
\hat{D}_{VA} &\sim \left( \Gamma - 1 \right) X_V, \;\;\;\;\;\;\;\;\;\;\;\;\;\, \hat{D}_{VV} \sim \left[\left(\Gamma - 1 \right)X_A + 1 \right].
\end{split}
\end{equation}
The boundary conditions are
\begin{equation}
\label{bc}
\frac{\partial X_A}{\partial x}\bigg|_{x=\pm l} = \frac{\partial X_V}{\partial x}\bigg|_{x=\pm l} = 0,
\end{equation}
and the initial conditions are
\begin{equation}
\begin{aligned}
\label{ic}
X_A(x,0) = \left\{
     \begin{array}{lcr}
      X_{A,\text{ini}} & \text{if} & -1 < \,\, x < 0,\\
      0 & \text{if} & 0 < \,\, x < 1,
     \end{array}
   \right. \;\;\;&\;\;\;
X_B(x,0) = \left\{
     \begin{array}{lcr}
      0 & \text{if} & -1 < \,\, x < 0,\\
      X_{B,\text{ini}} & \text{if} & 0 < \,\, x < 1,
     \end{array}
   \right. \\
   \; \\
   X_V(x,0) &= X_{V,\text{ini}},
\end{aligned}
\end{equation}
where $X_{A,\text{ini}}$, $X_{B,\text{ini}}$, and $X_{V,\text{ini}}$ denote the constant initial mole fractions of material A, B, and vacancies, respectively, and $X_{i,\text{ini}} = 1 - X_{V,\text{ini}}$, for $i=$ A, B.

In the following section we will show that the CA model reduces to \eqref{gov-eq-1-mathmodel}, \eqref{gov-eq-2-mathmodel} in the appropriate limit.

\section{Cellular automaton model}
\label{sec-CA}

In order to define the cellular automaton model for the one-dimensional problem discussed in \S \ref{sec-model-reduced} we are going to define a two-dimensional space of size $N\times N$ that is partitioned into two-dimensional $1\times 1$ cells. This would correspond to a two dimensional lattice of $N\times N$ atoms, in which each cell corresponds to one atomic site. Thus, the list of states for each cell in the CA model are ``A atom'', ``B atom'', and ``vacancy V''. The grid is considered to be periodic on the top/bottom edges. The model presented here will be asynchronous, that is, at each time step only one cell will be picked to apply the state change rule. Since physically atomic diffusion happens via vacancy exchange, it makes sense that the only cells in our CA grid that change state are those situated next to a vacancy cell and the vacancy cells themselves. For this reason, at each time step we only pick cells that represent vacancies to apply the change of state rule. Moreover, the choice of which vacancy cell is picked is done at random. We define the local neighbourhood of a cell as all the cells that surround it. Thus, each cell has eight neighbours, except the ones on the left and right columns on the grid, which can have five or three (corners) neighbours. We pick one of these neighbours at random and then apply the state change rule, which is defined as follows. If the neighbour cell picked is an A cell, we will proceed to exchange the states of the vacancy and A cell, and so A has moved. If the neighbour cell picked is a B cell, the probability of exchanging states with the vacancy cell is defined to be $1/\Gamma$. This will capture the physical feature in the model of B being $\Gamma$ times slower than A. Finally, if the neighbour cell picked is a vacancy, no change of states is applied.

Our interest now is to find the continuum limit of the asynchronous cellular automaton model we have described. Let us define the fraction of V cells, A cells, and B cells in the whole grid as
\begin{equation}
\bar{V} = \frac{N_V}{N^2}, \qquad \bar{A} = \frac{N_A}{N^2}, \qquad \bar{B} = \frac{N_B}{N^2},
\end{equation}
where $N_i$ is the number of cells of state $i$. Let us pick a square subgrid of size $\sqrt{N} \times \sqrt{N}$, and name it the $(i,j)$ subgrid (see Figure \ref{figure-subgrid}). Inside it, we define the following three functions,
\begin{itemize}
\setlength\itemsep{0em}
\item $\bar{V}\ensuremath{_{i,\,j}}$, the fraction of V cells in the $(i,j)$ subgrid;
\item $\bar{A}\ensuremath{_{i,\,j}}$, the fraction of A cells in the $(i,j)$ subgrid;
\item $\bar{B}\ensuremath{_{i,\,j}}$, the fraction of B cells in the $(i,j)$ subgrid.
\end{itemize}
All three functions above are dependent on space and time. Note that the choice of $M=\sqrt{N}$ of the subgrid is arbitrary. We only need $\lim_{N\rightarrow \infty} M/N = 0$ in a suitable manner.

We wish to study the evolution in one time step of the fraction number of V cells, A cells and B cells inside the $(i,j)$ subgrid. Let $\bar{A}\ensuremath{_{i,\,j}}^n$, $\bar{V}\ensuremath{_{i,\,j}}^n$, $\bar{B}\ensuremath{_{i,\,j}}^n$ be the fraction number of A cells, vacancies and B cells, respectively, in the $(i,j)$ subgrid at time step $n$. The aim is to compute $\bar{A}\ensuremath{_{i,\,j}}\ensuremath{^{n+1}}$, $\bar{V}\ensuremath{_{i,\,j}}\ensuremath{^{n+1}}$ and $\bar{B}\ensuremath{_{i,\,j}}\ensuremath{^{n+1}}$. In the next two sections we will discuss the change of A cells and V cells in one time step inside the  $(i,j)$ subgrid, respectively. We will omit the case of B cells since by conservation it can be found from A and V.

\begin{figure}[!htb]
\centering
\begin{tikzpicture}[scale=0.5, transform shape]
\draw (-8.5,-5) -- (8.5,-5);
\draw (-8.5,-3) -- (8.5,-3);
\draw (-8.5,-1) -- (8.5,-1);
\draw (-8.5,1) -- (8.5,1);
\draw (-8.5,5) -- (8.5,5);
\draw (-5,6.5) -- (-5,-8.5);
\draw (-3,6.5) -- (-3,-8.5);
\draw (-1,6.5) -- (-1,-8.5);
\draw (1,6.5) -- (1,-8.5);
\draw (3,6.5) -- (3,-8.5);
\draw (5,6.5) -- (5,-8.5);

\draw[line width=0.75mm] (-5.5,3) -- (5.5,3);
\draw[dashed, line width=0.75mm] (-5.5,3) -- (-6.5,3);
\draw[dashed, line width=0.75mm] (5.5,3) -- (6.5,3);
\draw[line width=0.75mm] (6.5,3) -- (7,3);
\draw[line width=0.75mm] (7,3) -- (7,-5.5);
\draw[dashed, line width=0.75mm] (7,-5.5) -- (7,-6.5);
\draw[line width=0.75mm] (7,-6.5) -- (7,-7);
\draw[line width=0.75mm] (7,-7) -- (6.5,-7);
\draw[dashed,line width=0.75mm] (6.5,-7) -- (5.5,-7);
\draw[line width=0.75mm] (5.5,-7) -- (-5.5,-7);
\draw[dashed, line width=0.75mm] (-5.5,-7) -- (-6.5,-7);
\draw[line width=0.75mm] (-6.5,-7) -- (-7,-7);
\draw[line width=0.75mm] (-7,-6.5) -- (-7,-7);
\draw[dashed,line width=0.75mm] (-7,-6.5) -- (-7,-5.5);
\draw[line width=0.75mm] (-7,3) -- (-7,-5.5);
\draw[line width=0.75mm] (-7,3) -- (-6.5,3);
\draw[line width=0.75mm] (-7,3) -- (-7,5.5);
\draw[dashed, line width=0.75mm] (-7,5.5) -- (-7,6.5);
\draw[line width=0.75mm] (7,3) -- (7,5.5);
\draw[dashed, line width=0.75mm] (7,5.5) -- (7,6.5);
\draw[line width=0.75mm] (-7,-7) -- (-7,-7.5);
\draw[dashed, line width=0.75mm] (-7,-7.5) -- (-7,-8.5);
\draw[line width=0.75mm] (7,-7) -- (7,-7.5);
\draw[dashed, line width=0.75mm] (7,-7.5) -- (7,-8.5);
\draw[line width=0.75mm] (-7,3) -- (-7.5,3);
\draw[dashed, line width=0.75mm] (-7.5,3) -- (-8.5,3);
\draw[line width=0.75mm] (7,3) -- (7.5,3);
\draw[dashed, line width=0.75mm] (7.5,3) -- (8.5,3);
\draw[line width=0.75mm] (-7,-7) -- (-7.5,-7);
\draw[dashed, line width=0.75mm] (-7.5,-7) -- (-8.5,-7);
\draw[line width=0.75mm] (7,-7) -- (7.5,-7);
\draw[dashed, line width=0.75mm] (7.5,-7) -- (8.5,-7);

\node at (0, 2) {\Huge V};
\node at (-4.75, -6) {$(i,j)$ subgrid};
\node at (-5.1, 5.4) {$(i,j+1)$ subgrid};
\node at (-5.1, -8) {$(i,j-1)$ subgrid};
\node at (-9, 2) {$(i-1,j)$ subgrid};
\node at (9, 2) {$(i+1,j)$ subgrid};

\node at (0, -9) {$\leftarrow \sqrt{N} \rightarrow$};
\node at (9, -1.5) {$\uparrow$};
\node at (9, -2) {$\sqrt{N}$};
\node at (9, -2.5) {$\downarrow$};

\draw [line width=0.75mm, red] (-3,5) -- (3,5);
\draw [line width=0.75mm, red] (-3,-1) -- (3,-1);
\draw [line width=0.75mm, red] (-3,5) -- (-3,-1);
\draw [line width=0.75mm, red] (3,5) -- (3,-1);

\end{tikzpicture}
\caption{Sketch of the subgrid set-up. In red the local neighbourhood of a V cell is shown.}
\label{figure-subgrid}
\end{figure}

\subsection{A cells}

There are two factors that can affect the amount of A cells in the $(i,j)$ subgrid in one time step: either a vacancy of the subgrid is able to exchange places with an A cell of a neighbouring subgrid (that adds an A cell), or a vacancy from a neighbouring subgrid is able to exchange places with an A cell in the $(i,j)$ subgrid (that removes an A cell). Thus, to add one A cell it is necessary that at the $n$-th time step
\begin{enumerate}[label={\{1.\arabic*\}},leftmargin=15mm]
\setlength\itemsep{0em}
\item a vacancy V inside the $(i,j)$ subgrid is picked;
\item said V is on one of the edges of the $(i,j)$ subgrid;
\item the cell picked to do the exchange is an A cell and is in one of the neighbouring $(i+1,j)$, $(i-1,j)$, $(i,j+1)$, $(i,j-1)$ subgrids.
\end{enumerate}
Similarly, to remove one A cell, it is necessary that at the $n$-th time step
\begin{enumerate}[label={\{2.\arabic*\}},leftmargin=15mm]
\setlength\itemsep{0em}
\item a vacancy V in one of the neighbouring $(i+1,j)$, $(i-1,j)$, $(i,j+1)$, $(i,j-1)$ subgrids is picked;
\item said V is on one of the edges of the subgrid it is on (the one neighbouring the $(i,j)$ subgrid);
\item the cell picked to do the exchange is an A cell and is in the $(i,j)$ subgrid.
\end{enumerate}

The probabilites of the events mentioned above are obtained via standard probability theory under the assumption that A and V are uniformly distributed in the subdomain,
\begin{align}
P(\{1.1\}) &= \frac{N \bar{V}\ensuremath{_{i,\,j}}}{N^2 \bar{V}}, \qquad P(\{1.2\}) = 4\frac{\sqrt{N}}{N}, \\
\label{A-exchange-prob}
P(\{1.3\}) &= \frac{1}{4}\left(\bar{A}\ensuremath{_{i+1,\,j}} + \bar{A}\ensuremath{_{i-1,\,j}} + \bar{A}\ensuremath{_{i,\,j+1}} + \bar{A}\ensuremath{_{i,\,j-1}} \right),\end{align}
and
\begin{align}
P(\{2.1\}) &= \frac{N}{N^2 \bar{V}}\left( \bar{V}\ensuremath{_{i-1,\,j}} + \bar{V}\ensuremath{_{i+1,\,j}} + \bar{V}\ensuremath{_{i,\,j-1}} + \bar{V}\ensuremath{_{i,\,j+1}} \right), \\
P(\{2.2\}) &= \frac{\sqrt{N}}{N}, \qquad P(\{2.3\}) = \bar{A}\ensuremath{_{i,\,j}}.
\end{align}

Now, the fraction of A cells at the next time step $\bar{A}\ensuremath{_{i,\,j}}\ensuremath{^{n+1}}$ is just the fraction of A cells at the current time step $\bar{A}\ensuremath{_{i,\,j}}^n$ plus the probability of adding an A cell into the subgrid, P(\{1.1\}) $\times$ P(\{1.2\}) $\times$ P(\{1.3\}), minus the probability of removing an A cell into the subgrid, P(\{2.1\}) $\times$ P(\{2.2\}) $\times$ P(\{2.3\}),
\begin{align}
\label{A-change}
\begin{split}
\bar{A}\ensuremath{_{i,\,j}}\ensuremath{^{n+1}} = &\bar{A}\ensuremath{_{i,\,j}}^n + \frac{1}{N} \left( \frac{N \bar{V}\ensuremath{_{i,\,j}}^n}{N^2 \bar{V}} \right) \left( 4\frac{\sqrt{N}}{N} \right) \frac{1}{4}\left(\bar{A}\ensuremath{_{i+1,\,j}}^n + \bar{A}\ensuremath{_{i-1,\,j}}^n + \bar{A}\ensuremath{_{i,\,j+1}}^n + \bar{A}\ensuremath{_{i,\,j-1}}^n \right)
\\
&-\frac{1}{N} \left( \frac{N}{N^2 \bar{V}}\left( \bar{V}\ensuremath{_{i-1,\,j}}^n + \bar{V}\ensuremath{_{i+1,\,j}}^n + \bar{V}\ensuremath{_{i,\,j-1}}^n + \bar{V}\ensuremath{_{i,\,j+1}}^n \right) \right) \left( \frac{\sqrt{N}}{N} \right) \bar{A}\ensuremath{_{i,\,j}}^n.
\end{split}
\end{align}

\subsection{V cells}

As in the previous case, there are two things can affect the amount of V cells in the $(i,j)$ subgrid: either an A or B cell in the $(i,j)$ subgrid is able to exchange places with a V cell in one of the neighbouring subgrids (that adds a V cell), or a vacancy from the $(i,j)$ subgrid is able to exchange places with an A or B cell in one of the neighbouring subgrids (that removes a V cell). To add one vacancy, it is necessary that at the $n$-th time step
\begin{enumerate}[label={\{3.\arabic*\}},leftmargin=15mm]
\setlength\itemsep{0em}
\item a vacancy V in one of the neighbouring $(i+1,j)$, $(i-1,j)$, $(i,j+1)$, $(i,j-1)$ subgrids is picked;
\item said V is on one of the edges of the subgrid it is on (the one neighbouring the $(i,j)$ subgrid);
\item the cell picked to do the exchange is on the $(i,j)$ subgrid and an exchange actually occurs. Recall that if the exchange is with an A cell the probability of movement is $1$ whereas if the exchange is with a B cell said probability is $1/\Gamma$.
\end{enumerate}
Similarly, to remove a vacancy, it is necessary that at the $n$-th time step
\begin{enumerate}[label={\{4.\arabic*\}},leftmargin=15mm]
\setlength\itemsep{0em}
\item a vacancy V inside the $(i,j)$ subgrid is picked;
\item said V is on one of the edges of the $(i,j)$ subgrid;
\item the cell picked to do the exchange is on one of the neighbouring $(i+1,j)$, $(i-1,j)$, $(i,j+1)$, $(i,j-1)$ subgrids and an exchange actually occurs.
\end{enumerate}
It is assumed that the vacancy concentration is low enough so that the probability of picking a vacancy to do the exchange is negligible. Consequently we only define the probabilities
\begin{align}
P(\{3.3\}) &=  1-\bar{B}_{i,j} \frac{\Gamma - 1}{\Gamma},\\
P(\{4.3\}) &= 1-\frac{1}{4} (\bar{B}\ensuremath{_{i+1,\,j}} + \bar{B}\ensuremath{_{i-1,\,j}} + \bar{B}\ensuremath{_{i,\,j+1}} + \bar{B}\ensuremath{_{i,\,j+1}}) \frac{\Gamma - 1}{\Gamma}.
\end{align}
This allows us to write
\begin{align}
\label{V-change}
\begin{split}
\bar{V}\ensuremath{_{i,\,j}}\ensuremath{^{n+1}} = &\bar{V}\ensuremath{_{i,\,j}} + \frac{1}{N} \left( \frac{N}{N^2 \bar{V}}\left( \bar{V}\ensuremath{_{i-1,\,j}}^n + \bar{V}\ensuremath{_{i+1,\,j}}^n + \bar{V}\ensuremath{_{i,\,j-1}}^n + \bar{V}\ensuremath{_{i,\,j+1}}^n \right) \right) \left( \frac{\sqrt{N}}{N} \right) \left( 1 - \bar{B}_{i,j}^n \frac{\Gamma - 1}{\Gamma} \right) \\
&-\frac{1}{N} \left( \frac{N V\ensuremath{_{i,\,j}}^n}{N^2 \bar{V}} \right) \left( 4 \frac{\sqrt{N}}{N}\right) \left( 1 - \frac{1}{4} (\bar{B}\ensuremath{_{i+1,\,j}}^n + \bar{B}\ensuremath{_{i-1,\,j}}^n + \bar{B}\ensuremath{_{i,\,j+1}}^n + \bar{B}\ensuremath{_{i,\,j+1}}^n)\frac{\Gamma - 1}{\Gamma} \right).
\end{split}
\end{align}

\subsection{Limit $N \rightarrow \infty$}
\label{sec-limit}

The objective of the present section is to determine whether in the limit $N \rightarrow \infty$ the CA model of the previous section reduces to the diffusion model in equations \eqref{gov-eq-1-mathmodel}-\eqref{gov-eq-2-mathmodel}.
Rearranging equation \eqref{A-change} leads to
\begin{align}
\label{A-change-limit}
\begin{split}
\bar{A}\ensuremath{_{i,\,j}}\ensuremath{^{n+1}} = &\bar{A}\ensuremath{_{i,\,j}}^n + \frac{\sqrt{N}}{N^2 \bar{V}} \left( \frac{\bar{V}\ensuremath{_{i,\,j}}^n \left( \bar{A}\ensuremath{_{i-1,\,j}}^n + \bar{A}\ensuremath{_{i+1,\,j}}^n + \bar{A}\ensuremath{_{i,\,j-1}}^n + \bar{A}\ensuremath{_{i,\,j+1}}^n \right)}{N} \right. \\
&\left.-\frac{\left( \bar{V}\ensuremath{_{i-1,\,j}}^n + \bar{V}\ensuremath{_{i+1,\,j}}^n + \bar{V}\ensuremath{_{i,\,j-1}}^n + \bar{V}\ensuremath{_{i,\,j+1}}^n \right) \bar{A}\ensuremath{_{i,\,j}}^n}{N} \right).
\end{split}
\end{align}
We may equate this to a standard finite difference form by choosing $\sqrt{N} = \Delta x$ and $\Delta \hat{t} = \frac{1}{N^{3/2} \bar{V}}$,
\begin{align}
\label{A-change-limit-2}
\begin{split}
\frac{\bar{A}\ensuremath{_{i,\,j}}\ensuremath{^{n+1}} - \bar{A}\ensuremath{_{i,\,j}}^n}{\Delta \hat{t}} = &\left( \frac{\bar{V}\ensuremath{_{i,\,j}}^n \left( \bar{A}\ensuremath{_{i-1,\,j}}^n + \bar{A}\ensuremath{_{i+1,\,j}}^n + \bar{A}\ensuremath{_{i,\,j-1}}^n + \bar{A}\ensuremath{_{i,\,j+1}}^n \right)}{\Delta x^2} \right. \\
&\left.-\frac{\left( \bar{V}\ensuremath{_{i-1,\,j}}^n + \bar{V}\ensuremath{_{i+1,\,j}}^n + \bar{V}\ensuremath{_{i,\,j-1}}^n + \bar{V}\ensuremath{_{i,\,j+1}}^n \right) \bar{A}\ensuremath{_{i,\,j}}^n}{\Delta x^2} \right).
\end{split}
\end{align}
Note that $\Delta \hat{t}$ is not a time, rather we are just conforming to finite difference notation.

Similarly, equation \eqref{V-change} may be expressed as
\begin{align}
\label{V-change-limit-2}
\begin{split}
\frac{\bar{V}\ensuremath{_{i,\,j}}\ensuremath{^{n+1}} - \bar{V}\ensuremath{_{i,\,j}}^n}{\Delta \hat{t}} = &\frac{\left( \bar{V}\ensuremath{_{i-1,\,j}}^n + \bar{V}\ensuremath{_{i+1,\,j}}^n + \bar{V}\ensuremath{_{i,\,j-1}}^n + \bar{V}\ensuremath{_{i,\,j+1}}^n \right)   \left( 1 - \bar{B}_{i,j}^n \frac{\Gamma - 1}{\Gamma} \right)}{\Delta x^2} \\
&-  \frac{\bar{V}\ensuremath{_{i,\,j}}^n \left( 4 - (\bar{B}\ensuremath{_{i+1,\,j}}^n + \bar{B}\ensuremath{_{i-1,\,j}}^n + \bar{B}\ensuremath{_{i,\,j+1}}^n + \bar{B}\ensuremath{_{i,\,j+1}}^n)\frac{\Gamma - 1}{\Gamma} \right)}{\Delta x^2}.
\end{split}
\end{align}

Taking the limit $N \rightarrow \infty$ in equations \eqref{A-change-limit-2} and \eqref{V-change-limit-2}, changing the notation to $A\ensuremath{_{i,\,j}} \equiv X_A$, $B\ensuremath{_{i,\,j}} \equiv X_B$, $V\ensuremath{_{i,\,j}} \equiv X_V$, and substituting $X_B = 1 - X_A - X_V$ gives
\begin{align}
\label{A-change-limit-4}
\frac{\partial X_A}{\partial \hat{t}} &= X_V \nabla^2 X_A - X_A \nabla^2 X_V,\\
\label{V-change-limit-4}
\Gamma \frac{\partial X_V}{\partial \hat{t}} &=  -(\Gamma - 1) X_V \nabla^2 X_A + \left( 1 + (\Gamma - 1)X_A \right)\nabla^2 X_V.
\end{align}

Finally, setting $t = \hat{t}/\Gamma$, we recover equations \eqref{gov-eq-1-mathmodel}-\eqref{gov-eq-2-mathmodel} and the CA model does indeed reduce to the continuous diffusion model.

\section{Results}

In this section we present the results of the CA model. For the simulations we pick a square grid $N \times N$, where $N=200$. At the first time step, the first 100 columns are A cells, and the remaining 100 columns are B cells. Then we randomly distribute 2000 vacancies (equivalent to 5\% of the total number of cells) throughout the whole grid. A simulation is then run for $1.85 \times 10^9$ steps which is sufficient to allow for significant change in the distribution of cells.

Figures \ref{figure-CA-mesh-different-times}(a), \ref{figure-CA-mesh-different-times}(c), \ref{figure-CA-mesh-different-times}(e) show the distribution of material and vacancies when $\Gamma = 1.5$, that is, the diffusion rates between A and B are similar. Throughout the process vacancies are well distributed in the domain and, by the time $t=t_3$, the system appears close to equilibrium. Figures \ref{figure-CA-mesh-different-times}(b), \ref{figure-CA-mesh-different-times}(d), \ref{figure-CA-mesh-different-times}(f) show the corresponding evolution when A diffuses much faster than B. In Figure \ref{figure-CA-mesh-different-times}(b) we see a greater motion of A to the right than the one observed in Figure \ref{figure-CA-mesh-different-times}(a). However this also means that vacancies accumulate on the left. In Figure \ref{figure-CA-mesh-different-times}(d) it is clear that the vacancy concentration on the right is low, which acts to slow down the diffusion. This is clear from the final figure, Figure \ref{figure-CA-mesh-different-times}(f), which is far from equilibrium. This may seem counter-intuitive; A diffuses much faster here than in Figure \ref{figure-CA-mesh-different-times}(e) but it clearly ends up moving slower. This is a result of the initial rapid movement of A, bringing a high proportion of vacancies to the left and so hindering further movement. A similar result was noticed in the continuum model of \cite{Ribera2017}, where a very fast diffuser ends up redistributing more slowly.

\begin{figure}[!htb]
\begin{center}
\begin{subfigure}[b]{0.49\textwidth}
\begin{center}
\includegraphics[width=0.78\textwidth]{CA_G1p5_grid_1-eps-converted-to.pdf}
\end{center}
\vskip-4mm
\caption{$\Gamma = 1.5$, $t = t_1$.}
\end{subfigure}
\begin{subfigure}[b]{0.49\textwidth}
\begin{center}
\includegraphics[width=0.78\textwidth]{CA_G10_grid_1-eps-converted-to.pdf}
\end{center}
\vskip-4mm
\caption{$\Gamma = 10$, $t = t_1$.}
\end{subfigure}

\vskip5mm

\begin{subfigure}[b]{0.49\textwidth}
\begin{center}
\includegraphics[width=0.78\textwidth]{CA_G1p5_grid_2-eps-converted-to.pdf}
\end{center}
\vskip-4mm
\caption{$\Gamma = 1.5$, $t = t_2$.}
\end{subfigure}
\begin{subfigure}[b]{0.49\textwidth}
\begin{center}
\includegraphics[width=0.78\textwidth]{CA_G10_grid_2-eps-converted-to.pdf}
\end{center}
\vskip-4mm
\caption{$\Gamma = 10$, $t = t_2$.}
\end{subfigure}

\vskip5mm

\begin{subfigure}[b]{0.49\textwidth}
\begin{center}
\includegraphics[width=0.78\textwidth]{CA_G1p5_grid_3-eps-converted-to.pdf}
\end{center}
\vskip-4mm
\caption{$\Gamma = 1.5$, $t = t_3$.}
\end{subfigure}
\begin{subfigure}[b]{0.49\textwidth}
\begin{center}
\includegraphics[width=0.78\textwidth]{CA_G10_grid_3-eps-converted-to.pdf}
\end{center}
\vskip-4mm
\caption{$\Gamma = 10$, $t = t_3$.}
\end{subfigure}

\caption{Resulting grid $200 \times 200$ for different times obtained with the simulation of the CA model. Red denotes A atom cells. Blue denotes B atom cells. Yellow denotes vacancy cells. $t_1 < t_2 < t_3$.}
\label{figure-CA-mesh-different-times}
\end{center}
\end{figure}

In Figure \ref{figure-CA-vs-FullNumerics} we compare the concentration of vacancies given by the CA model and that given by the continuum model in equations \eqref{gov-eq-1-mathmodel}-\eqref{gov-eq-2-mathmodel} at different times for two values of ratio between the diffusion rates; $\Gamma = 1.5$ (Figure \ref{figure-CA-vs-FullNumerics}(a) and \ref{figure-CA-vs-FullNumerics}(b)) and $\Gamma = 10$ (Figure \ref{figure-CA-vs-FullNumerics}(c) and \ref{figure-CA-vs-FullNumerics}(d)). The numerical solution of \eqref{gov-eq-1-mathmodel}-\eqref{gov-eq-2-mathmodel} is standard and defined in \cite{Ribera2017}. The CA results come from an average of simulations. To achieve this we define
\begin{equation}
\label{average-V}
V_j^n = \frac{\bar{V}}{I} \sum_{i=1}^I \left\{N_V^i \right\}^n_j
\end{equation}
where $\left\{N_V^i \right\}^n_j$ is the number of vacancies at the $j$-th column on the $n$-th time step, the superscript $i$ is used to distinguish different simulations, and $I$ is the total number of CA simulations. The variable $V_j^n$ denotes the average concentration of vacancies at the $j$-th column on the $n$-th time step. The results shown in Figure \ref{figure-CA-vs-FullNumerics} correspond to $I=10$.

To be able to compare the variable $V_j^n$ and the numerical solution $X_V$ of \eqref{gov-eq-1-mathmodel}-\eqref{gov-eq-2-mathmodel} we need to find a correspondence between the discrete time in the CA simulation and the continuous time in the PDE system. Let $n$ be the time step that needs to be transformed to a continuum time $t_n$. Then
\begin{equation}
t_n = n \frac{\Delta \hat{t}}{\Gamma} t_s^{-1},
\end{equation}
where $t_s$ is the time scale defined in \cite{Ribera2017}. It is defined as $t_s =  \lambda a^2 \Gamma_B/(\epsilon l^2)$,
where $\epsilon = 2 X_{V,\text{ini}}$, $\lambda$ is a geometric factor, $a$ is the lattice constant and $\Gamma_B$ is the hop frequency of species B. Finally, $l$ is the length of the one-dimensional bar and it is found via the lattice site density $\rho$ and the number of cells in our CA grid,
\begin{equation}
l = \left(\frac{N^3}{\rho}\right)^{1/3},
\end{equation}
which gives $l=51$ nm for the case presented in Figure \ref{figure-CA-vs-FullNumerics}. Using the parameters values shown in \cite[Table 1]{Ribera2017}, we find that $t_s \approx 1.05 \times 10^{8}$. There is only one issue remaining to be treated before being able to compare the two solutions. When solving the continuum model, the initial condition for vacancies is given by \cite{Ribera2017}
\begin{equation}
X_{V,1}(x,0) = \frac{1+(\Gamma-1)X_A(x,0)}{2+(\Gamma-1)X_{A,\text{ini}}},
\end{equation}
where $X_{V,1} = X_V/\epsilon$. We do this rescaling to be able to keep track of the evolution of vacancies since this number is usually very small compared to the concentration of species A and B; numerically, this may cause problems. The initial condition for the CA model corresponds to $X_{V,1}(x,0) \approx 0.5$. This merely means that if $n_0$ is the initial time step in the CA model, the actual $n_{t_0}$ that corresponds to the initial time $t_0$ is $n_{t_0} = r n_0$, where $r > 1$ (see Figure \ref{sketch-conversion}). To find this rescaling factor $r$ we minimise the least-squares error between the continuum model and the discrete set of data $V_j^{n_f}$ at $n_f=1.85\times 10^9$. For $\Gamma = 10$, we find that $r = 26.171$; for $\Gamma = 1.5$, $r = 21.528$.
The comparison of $V_j^{n}$, for $j=1,\dots,\,N$ and $X_{V,1}(t,x)$, with $x\in [-1,\,1]$, is now well defined.
\begin{wrapfigure}{R}{7.5cm}
\centering
\begin{tikzpicture}
\node at (0,0.6) {CA discrete time};
\node at (4,0.6) {Continuous time};
\node at (0,0) {$n_0$};
\draw[arrows=->] (0.4,0) -- (3.7,0);
\draw[arrows=->] (0,-0.2) -- (0,-1.7);
\node at (4,0) {$\nexists$};

\node at (0,-2) {$n_{t_0}$};
\draw[arrows=->] (0.4,-2) -- (3.7,-2);
\node at (4,-2) {$t_0$};

\draw[arrows=->,dashed] (0.4,-0.2) -- (3.7,-1.8);

\node at (2.9,-0.85) {$r n_0 = n_{t_0}$};

\end{tikzpicture}
\caption{Time conversion sketch.}
\label{sketch-conversion}
\label{sketch-conversion}
\end{wrapfigure}

In Figure \ref{figure-CA-vs-FullNumerics} we compare results for the vacancy concentration from the continuum model of equations \eqref{gov-eq-1-mathmodel}-\eqref{gov-eq-2-mathmodel} and the average result of 10 simulations via equation \eqref{average-V}. First, we note what was observed in the previous figures, when $\Gamma = 1.5$ the vacancy concentration is relatively constant. When $\Gamma = 10$ vacancies concentrate on the left, thus slowing the movement of the fast diffuser. All figures show good agreement, even when $\Gamma = 1.5$. The most noticeable discrepancies occur at small times, near the ends $x=\pm 1$, where the continuum model indicates greater movement from the initial condition ($X_V = 0.5$). This is not surprising, continuum diffusion models typically allow for motion throughout the domain even though in reality extreme points may not be feeling any effect. Hence we expect that for small times the CA model is more realistic, and the agreement improves with time.

\begin{figure}[!htb]
\begin{center}
\begin{subfigure}[b]{0.49\textwidth}
\includegraphics[width=\textwidth=1]{CA_vs_Num_G1p5_t1-eps-converted-to.pdf}
\caption{$\Gamma = 1.5$, $t = 0.0378$.}
\end{subfigure}
\begin{subfigure}[b]{0.49\textwidth}
\includegraphics[width=\textwidth=1]{CA_vs_Num_G1p5_t2-eps-converted-to.pdf}
\caption{$\Gamma = 1.5$, $t = 0.2519$.}
\end{subfigure}
\begin{subfigure}[b]{0.49\textwidth}
\includegraphics[width=\textwidth=1]{CA_vs_Num_G10_t1-eps-converted-to.pdf}
\caption{$\Gamma = 10$, $t = 0.0574$.}
\end{subfigure}
\begin{subfigure}[b]{0.49\textwidth}
\includegraphics[width=\textwidth=1]{CA_vs_Num_G10_t2-eps-converted-to.pdf}
\caption{$\Gamma = 10$, $t = 0.3263$.}
\end{subfigure}
\caption{Dotted line is obtained by joining the discrete normalised vacancy average concentration $V_j^n$ obtained via the CA model (equation \eqref{average-V}). Solid line represents the numerical solution of the continuum model described by equations \eqref{gov-eq-1-mathmodel}-\eqref{gov-eq-2-mathmodel}.}
\label{figure-CA-vs-FullNumerics}
\end{center}
\end{figure}

\section{Conclusions}

The goal of this paper was to develop a cellular automaton model to describe binary diffusion in solids. To do this we chose to employ an asynchronous model and used basic rule definitions to update the cell states, based on the physics that drive the process.

In the limit of large number of cells we showed that this CA model reduces to a particular form of the continuum model developed in \cite{Ribera2017}. This was verified in the Results section by comparing the CA and the continuum model. This opens the possibility of designing various scenarios in a very simple fashion and then just take the limit to obtain the continuum model to do a more accurate analysis.

An interesting result to come out of this work is that when the system has a very fast diffuser it can lead, overall, to a slower diffusion process. This occurs because the initial fast diffusion acts to move nearly all vacancies to one side, so restricting further vacancy exchange and so movement. Also, at small times the nature of the continuum model permitted movement throughout the domain when in practice this may not occur. The CA model showed less movement near the extreme points at small times. This seems more physically realistic, hence in this case the CA model may be preferable.

Obviously the CA model becomes increasingly cumbersome as the number of cells increases. For sufficiently large numbers a continuum model is clearly preferable. However, when the number of cells is small, for example when modelling nanoscale diffusion, CA models provide a powerful tool which may be more accurate than the continuum models.

\section*{Acknowledgements}

HR and TM acknowledge that the research leading to these results has received funding from ``la Caixa'' Foundation and has been partially funded by the CERCA Programme of the Generalitat de Catalunya. TM acknowledges Ministerio de Ciencia e Innovación grant MTM2014-5621. BW acknowledges an NSERC Canada research grant.

\clearpage