\section{Introduction}
Big data and machine learning are monumental to the furthering of science. With many existing large image data sets \citep{imagenet, mnist}, machine learning and deep learning have become extremely effective at classifying such data \citep{dai2021coatnet, an2020ensemble}. Data sets for classification consist of different object types, or classes, and the true class labels for the images are known as the "ground truth". The task of classifying images given the ground truth of a training data set is known as supervised classification. In this study we perform classification of greyscale images from the Caltech Spheromak Experiment \citep{bellan}, where plasma jets and spheromaks are created and photographed in a vacuum chamber. Initial classification of data from this experiment is desirable in that it provides a starting point for further investigation into machine recognition of the underlying physics in the experiment from optical images, as well as extrapolation from laboratory plasma jets to our understanding of astrophysical plasma jets.

Convolutional neural networks (CNNs) are known to excel at supervised image classification \citep[see][]{cnn_rev}, where they achieve state of the art accuracy on benchmarks such as the Imagenet benchmark \citep{imagenet_bm}. Feature engineering is another common practice in machine learning and classification, where hand-defined, or engineered, features of images are often used in classification tasks to differentiate between images. The values of these features are organized into an array, which is used as a feature vector representing an image in the feature space. This study classifies both engineered feature vectors and output feature vectors from CNNs through the angle between feature vectors or similarly, the cosine of the angle between feature vectors known as the cosine similarity. Figure 1 displays a mock example of how this is done.

\begin{figure}
  \centering
  \includegraphics[scale = 0.8]{fig1.png
  \caption{We conduct classification using the cosine similarity between feature vectors that are extracted from images. In order for classification to be conducted, the vectors must greatly align for images in the same class compared to images in different classes. Classification reduces to probing if a vector lies within a particular range in the binary classification case (shown above) or choosing the closest target region in the multi-class classification case. In this example, the similar image would be classified in the same class as the target image, while the dissimilar image would be classified out of the class.}
\label{fig:1}
\end{figure}

Cosine similarity has been used in face verification \citep{Nguyen2010CosineSM}, it can replace the dot product between output layers and weight vectors before activation layers in multi-layer networks \citep{luo2017cosine}, and it has been used to select robust feature subsets in genomic data analysis \citep{cs_gen}. Cosine similarity has also been used to adaptively refine clustering results in a pair-wise binary classification scheme \citep{dac}, and it is used in natural language processing to compare word embeddings \citep[see][]{young2018recent}. Cosine similarity has not yet been used to classify laboratory plasma images, and an algorithm that performs well using it to classify feature vectors is desirable given its simplicity. In this work we describe laboratory plasma image classification using cosine similarity.

In addition, this study uses AlexNet \citep{alexnet}, one of the most well known CNNs, to achieve notably accurate transfer learning. Furthermore, we use cosine similarity of image feature vectors as a loss function for our instance of AlexNet, to select relevant features, and to classify images directly. If cosine similarity shows reasonable results in the ability to characterize similarity/difference between vectors in different classes based on specific features, then it should also be able to perform in a simple classification algorithm using vector direction comparison akin to the mock example from Figure 1. This work validates this idea by selecting features that distinguish the feature vector directions for images in a specific class from others. We train an instance of AlexNet to increase the cosine similarity of output feature vectors of images in the same class while attempting to orthogonalize feature vectors of images in different classes. We find that using a simple direction comparison algorithm with the output vectors from this instance of AlexNet achieves classification accuracies notably higher than traditional transfer learning with cross entropy.

The remainder of this work is organized as: Methods and Results, where we give the details of the algorithms and image processing procedures while noting results obtained in using cosine similarity for feature selection and classification with AlexNet, Discussion, where the results and implications are assessed and future avenues are discussed, and Summary, where the key ideas from the paper are summarized.

\section{Methods and Results}
\subsection{Data Sets}
Our goal is to classify plasma jet images from the Caltech Spheromak Experiment \citep{scott}, where laboratory plasma jets are created and rich phenomena including kink instabilities are observed. Thus, our data sets consist of images sorted from this experiment. We create two data sets, one small data set (about 1000 images) with ground truth labels and one large data set (about 45000 images) with no ground truth labels.

For our small data set we sort the data into five classes, pertaining to five types of images that appear frequently in the experiment. Figure 2 shows a representative image from each class as an example, with a false colour map on each greyscale image. The classes are as follows: Haze, which corresponds to a general diffuse and spread-out plasma configuration, Spider, which corresponds to an initial eight plasma arcs tracing a background poloidal magnetic field occurring near the inception of a plasma jet, Column, a straight plasma jet showing little kink instability, Kink, an unstable column showing kink instability akin to a corkscrew shape, and Sphere, corresponding to anomalous spherical arrangements of pixels and plasma columns breaking off of the jet to form spheromaks. The number of images in each class for this small data set is 226, 111, 130, 295, and 390 for Haze, Spider, Column, Kink, and Sphere, respectively. We split each class in this set with a randomly distributed fractional split of 0.5, 0.25, and 0.25 and compile the corresponding splits of each class for training, validation, and test sets, respectively.

For our large data set, we separate images from the total data of about 300000 raw images using an instance of ResNet152 \citep{resnet} to distinguish between images based on quality.
We select 45000 high quality images for our large, unlabeled data set.
We set aside this large data set without labels for use in section 2.6 at the end of this work.
 For a thorough introduction to the experiment yielding these images and the theory involved, see \citep{bellan}.
\begin{figure}
  \centering
  \includegraphics[scale = 0.9]{fig2.png
  \caption{An example image from each of the five classes in the small data set. The classes are Haze, Spider, Column, Kink, and Sphere. Each class contains 226, 111, 130, 295, and 390 images respectively. Note that the images feature a false colour map.}
\label{fig:2}
\end{figure}

\subsection{Feature Selection Using Cosine Similarity}
To classify our small data set, we first turn to statistical and information theoretic features of images. We conduct feature extraction on our data set, inspired by the ability to view a greyscale image as a two dimensional intensity distribution. We transform each image into four statistical distributions, $P(x,y)$, $P(x)$, $P(y)$, and $P(v)$, for the collection of features from these distributions.  An example of these distributions obtained from an image can be seen in Figure 3, which displays an image, projections to the three distributions, $P(x,y)$, $P(x)$, and $P(y)$, and an example of the $P(v)$ distribution retrieved from the image. The quantitative description of each distribution is expressed in Appendix A. In order to obtain $P(x,y)$, the greyscale image matrix is normalized by dividing by the full intensity value of the image. $P(x)$ and $P(y)$ are obtained by summing along the y-axis and x-axis of the image, respectively, then normalizing by the total intensity of the image. $P(v)$ is obtained by taking a histogram of pixel values of an image.
\begin{figure}
  \centering
  \includegraphics[scale = 0.9]{fig3.png
  \caption{An example of an image from the kink class with the four distributions yielded to extract features. Note the image features a false colour map. $P(x,y)$, $P(x)$, and $P(y)$ are distributions of the intensity (left), while $P(v)$ is a histogram of pixels with particular values (right). Features are extracted from statistical and information theoretic quantities of these four distributions.}
\label{fig:3}
\end{figure}

With four distributions, we extract nineteen basic features from statistical and information theoretic quantities of the distributions, and we also include the total intensity value of each image as a feature. The quantitative description of each feature is expressed in Appendix A. The features are quite natural and easily calculable from an image and its distributions. They are as follows: the total intensity of the image (one feature), the raw mean and central variance of the $P(x)$, $P(y)$, and $P(v)$ distributions (six features), the third and fourth standardized moments of the $P(x)$, $P(y)$, and $P(v)$ distributions (six features), the entropy of each distribution (four features), the conditional entropies of the $P(x,y)$ distribution (two features), and the mutual information between the $P(x)$ and $P(y)$ distributions (one feature). These features can be computed on different pixel value ranges, allowing more features to be computed. We opt to compute the features on pixels valued from 200-300.

In order to search for features that can distinguish classes, we create "cosine similarity matrices" which display the pair-wise average cosine similarity of classes. Generally, for each pair of two classes, we compute the cosine similarity ($c.s. = \hat{u}\cdot\hat{v}$) for each pair of vectors between the classes and calculate the average. We then display the results in a cosine similarity matrix, where the cross sections display the average cosine similarity value obtained between classes. It is expected that a proper feature set will have a matrix with diagonals, which include the similarity of each class with itself, with higher values compared to the off-diagonals, such that a matrix with diagonals of one and off-diagonals of negative one would be the perfect result.

Specifically, we compute cosine similarity matrices yielded from all subsets of two features from the original twenty features and exhaustively search for a well behaved cosine similarity matrix indicating features geared towards a specific classification task. In Figure 4 we show a selection of a cosine similarity matrix from this search, with the mutual information of the $P(x)$ and $P(y)$ distributions, $I(x,y)$, and the mean of the $P(y)$ distribution for the features. From the cosine similarity of class one (spider) with itself compared to the cosine similarity of class one with the other four classes, these two features show potential for a binary, class one vs. non class one (spider vs non-spider) classification simply using comparison of feature vectors to a representative sample of class one vectors. Unfortunately, we do not see any matrices optimized for a full five-way classification of the data set from cosine similarity matrices of all possible feature vectors of length two, three, and four from our list of twenty features. Thus we conclude that these features are not good choices for our five-way plasma image classification goal.
\begin{figure}
  \centering
  \includegraphics[scale = 0.5]{fig4.png
  \caption{Example of a cosine similarity matrix indicating a potential feature selection for the task of class one or non class one (spider versus non-spider) binary classification. The two features used as components for vectors to construct the matrix were the mutual information, $I(x,y)$, and the mean of the $P(y)$ distribution. Classes are labeled as 0 through 4 corresponding to the classes shown in Figure 2. Since the matrix is normally symmetric, the bottom left half is set to the minimum value and crossed out.}
\label{fig:4}
\end{figure}
\subsection{A Binary Classification Algorithm Using Cosine Similarity of Feature Vectors}
We opt to construct a simple classifier using feature vector direction comparison for the following two reasons: one, to demonstrate the uses of cosine similarity for feature extraction by using the algorithm on the feature vectors that yielded Figure 4, and two, for the potential of the algorithm to be used with vectors that yield a more well behaved cosine similarity matrix and potential for multi-class classification.

Our binary classification algorithm is diagrammed in Figure 5 and outlined in the Algorithm 1 block. It is essentially a vector direction comparison algorithm. We first "train" the algorithm by selecting a representative sample of training images from the class. We extract feature vectors from each image in the training sample, and we compute the average cosine similarity of a large number of random pairs of feature vectors within the sample. Here the algorithm "learns" the similarity of images belonging to that class. The cosine similarity values calculated from "training" have a distribution with a training mean, $\mu$, and training standard deviation, $\sigma$. A test image, to be classified as in the class or not, is introduced to the algorithm, where its features are extracted and its cosine similarity is calculated with a satisfactory number (usually all) of the training samples. These new cosine similarity values have a test mean, $\mu^{*}$, which can be compared to the training mean. The training mean is representative of the average similarity an image in the class should have with its members, and the test mean represents the average similarity of the test image with the class. Classification is reduced to a threshold function, $\tau(\mu,\mu^{*},\sigma)$, which simply decides an image is in the class if $\mu^{*}$, the test mean, is greater than  $\mu - \gamma\sigma$. $\gamma$ is the only parameter featured in the algorithm.

\begin{algorithm}
\caption{: \textbf{Binary Classification With Cosine Similarity} \\
Given: iterations $N$, a number $\gamma$, a training set $U$, a test set $P$, and a vector-valued function $f(x)$}\label{alg:1}
\begin{algorithmic}
\State \:\:\:\:\:\:\:\:1.\:\:\:\:\textbf{for} $i$ in $N$: \Comment{"Training"}
\State \:\:\:\:\:\:\:\:2.\:\:\:\:\:\:\:\:\:\:\:\:$u \overset{\mathrm{iid}}{\sim}$ Unif$(U)$ ;  $v \overset{\mathrm{iid}}{\sim}$ Unif$(U)$
\State \:\:\:\:\:\:\:\:3.\:\:\:\:\:\:\:\:\:\:\:\:$x_i \xleftarrow{} \frac{f(u)\cdot f(v)}{\lVert f(u) \rVert \lVert f(v) \rVert}$
\State \:\:\:\:\:\:\:\:4.\:\:\:\:$\mu = \frac{1}{N}\sum\limits_{i=1}^N x_i$ ; $\sigma = \sqrt{\frac{\sum^{N}_{i=1} (x_i-\mu)^2}{N}}$
\State
\State \:\:\:\:\:\:\:\:5.\:\:\:\:\textbf{for} $p$ in $P$: \Comment{"Testing"}
\State \:\:\:\:\:\:\:\:6.\:\:\:\:\:\:\:\:\:\:\:\:\textbf{for} $u$ in $U$:
\State \:\:\:\:\:\:\:\:7.\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:$x_u \xleftarrow{} \frac{f(p)\cdot f(u)}{\lVert f(p) \rVert \lVert f(u) \rVert}$
\State \:\:\:\:\:\:\:\:8.\:\:\:\:\:\:\:\:\:\:\:\:$\mu^* = \frac{1}{\lVert U \rVert}\sum\limits_{u=1}^{\lVert U \rVert} x_u$
\State \:\:\:\:\:\:\:\:9.\:\:\:\:\:\:\:\:\:\:\:\:\textbf{if} $\mu^* > \mu - \gamma \sigma$:
\State \:\:\:\:\:\:\:\:10.\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:$p$ is in class
\State \:\:\:\:\:\:\:\:11.\:\:\:\:\:\:\:\:\:\:\textbf{else}:
\State \:\:\:\:\:\:\:\:12.\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:$p$ is NOT in class
\end{algorithmic}
\end{algorithm}

\begin{figure}
  \centering
  \includegraphics[scale = 0.25]{fig5.png
  \caption{A diagram of the binary classification algorithm based on cosine similarity. We start with a "training" stage, such that enough vectors are sampled from the spider class an approximate average cosine similarity, $\mu$, of members in the class is retrieved, along with a standard deviation, $\sigma$, of the distribution of cosine similarity values in the class. With this information we introduce a test image to the algorithm and compute its average cosine similarity, $\mu^{*}$, with the training samples. Then a threshold function, $\tau (\mu, \mu^{*}, \sigma)$, which considers an image a member in the class if $\mu^*$ is greater than  $\mu - \gamma\sigma$, determines whether the image is in the spider class or not. Note the plasma images feature a false colour map. This algorithm is used to achieve notable results (shown in Figure 6) with the features from Figure 4.}
\label{fig:5}
\end{figure}

This binary classification method is used on the small data set with the same features from Figure 4 as a feature vector of length two for each image. We implement the method with fifty training samples from the spider class and fifty test samples from each class totalling 250 images to be classified. We choose to sample 1000 random pairs in the training set for the computation of $\mu$, and every vector in the training set is compared with the test images for the computation of $\mu^*$. $\gamma$ is chosen to be $0.01$ for optimal results. The results of spider versus non-spider binary classification using these presets are shown through an extended confusion matrix in Figure 6, where the algorithm performs surprisingly well despite its simplicity. Given a spider image, the algorithm classifies it as a spider with 80\% accuracy. For the remaining classes, the algorithm correctly classifies images as non-spider with 87.5\% average accuracy. In order to test the need for large amounts of training data, the same confusion matrix is calculated using 100 training pairs instead of 1000 with which the algorithm correctly classifies spiders and non-spiders with 94\% and 83\% accuracy, respectively.
\begin{figure}
  \centering
  \includegraphics[scale = 0.9]{fig6.png
  \caption{Confusion matrix of the results using the binary classification algorithm on the features from Figure 4 for spider versus non-spider classification. The matrix displays fractional accuracy for each class in the first row and the fractional inaccuracy for each class presented in the second row. The algorithm achieves an average accuracy of 86\%, and the total number of images classified is 250.}
\label{fig:6}
\end{figure}

\subsection{Transfer Learning Using Cosine Embedding and Cross Entropy Loss}
The results of the binary classification algorithm motivated us to conduct transfer learning on our plasma data set using both the cosine embedding loss function and the cross entropy loss function. This was done to attempt the full five-way classification of our data set and to probe how the training on our data set with different loss functions, specifically one related to cosine similarity, would allow for a well behaved cosine similarity matrix from output vectors of the network.

We start AlexNet with the pretrained weights from pytorch \citep{pytorch}, which pertain to weights obtained upon the training of AlexNet on the Imagnet dataset \citep{imagenet}. The training, validation, and testing fractional split of our data is 0.5, 0.25, and 0.25. We use Amazon Web Service's Sagemaker for training along with hyperparameter optimization. We choose to optimize the learning rate (0.001 to 0.2), epoch number (5 to 100), and batch size (10 to 60) for cross entropy training and achieve five-way classification accuracies of 88\% for validation and 84\% for testing. For cosine similarity trained AlexNet, we choose the cosine embedding loss function, which works to maximize cosine similarity for objects in the same class and orthogonalizing objects in different classes. We implement the cosine embedding loss by using a hyperparameter, $\lambda$. In our classification, image labels are conventionally one-hot vectors, such as (0 0 1 0 0) corresponding to an image in class three from one to five. We then take predictions of the same form, where we simply take the maximum input in a vector, set it to one, and set the remaining inputs to zero. However, the cosine embedding loss between one-hot vectors results in gradients of one when the prediction is correct and zero when the prediction is incorrect. By setting all of the zeros in one-hot vectors to equal $\lambda$, we achieve the ability to conduct gradient descent, and thus, minimize the loss while also adding a parameter that tunes the direction of output vectors. If we set $\lambda$ to a constant negative one, we can perform gradient descent, but we heavily weight the direction of vectors through the inputs of negative one. We implement $\lambda$ with a range of  -0.4 to -0.05 while using the same hyperparameter ranges from cross entropy training for all other hyperparameters, and we achieve five-way classification accuracies of 89\% for validation and 82\% for testing.

For a visualization of the activation differences between the trained versions of AlexNet, we employ the Gradient-Weighted Class Activation Mapping method \citep{gradcam}. The method highlights key areas of an image that allow the network to recognize a class by using the incoming gradients to the last convolutional layer of a network. Since the last convolutional layer contains the last instance of spatial information before being flattened in the fully connected layers, it ideally contains the most relevant class-specific feature activations which contribute to the fully connected layers and classification.
The method is implemented using readily available code\footnotemark[1]\footnotetext[1]{https://github.com/utkuozbulak/pytorch-cnn-visualizations}, with results shown in Figure 7. We note that the cross entropy trained model prioritized less of the total shapes for recognition of the column and kink images opposed to the cosine similarity trained model, which uses the majority of the kink and column shapes for their recognition.
\begin{figure}
  \centering
  \includegraphics[scale = 1.5]{fig7.png
  \caption{Differences in the activations between the cosine similarity trained AlexNet (labeled as CS) versus the cross entropy trained AlexNet (CE). The use of Gradient-weighted Class Activation Mapping (Grad-CAM) was performed on the five representative plasma images from Figure 2 for both versions of AlexNet. GT stands for "Ground Truth". The second and third rows display the outputs of Grad-CAM method on the first row for the corresponding models. Note the ground truth images are displayed using a false colour map.}
\label{fig:7}
\end{figure}

In addition to training using the cross entropy and cosine embedding losses, we also perform gradient descent using the cosine embedding loss on the final 4096 dimensional ReLU layer of AlexNet, with the aim of yielding a model that retrieves outputs with a well-behaved cosine similarity matrix of embedded vectors for use with our classification algorithm. The general training process is shown in Figure 8.
\begin{figure}
  \centering
  \includegraphics[scale = 0.35]{fig8.png
  \caption{A diagram of the cosine embedding training on the final ReLU layer of AlexNet.}
\label{fig:8}
\end{figure}
We perform gradient descent with no classification results due to the final layer of the network no longer being a fully connected layer with outputs corresponding to classification.

We then use the three trained models and compute cosine similarity matrices on output vectors from the final ReLU layer to display the differences in the cosine similarity of the final vectors before classification. Note that we perform this using the images from the test set only. Figure 9 shows these results, which demonstrate that conventional transfer learning of AlexNet does not necessarily include an optimization of cosine similarity between the final ReLU layer for images in the same class. Training by performing gradient descent on the embeddings themselves potentially gives the model the ability to perform with our vector direction comparison algorithm, so much so that vector direction comparison using cosine similarity could potentially be used to conduct full five-way classification of the test data set.
\begin{figure}
  \centering
  \includegraphics[scale = 0.20]{fig9.png
  \caption{Cosine similarity matrices computed after performing transfer learning on AlexNet using a cosine embedding loss or cross entropy loss. Note the vectors used to construct the matrices were from the final ReLU layer of each network and that the matrices were computed using images from the test set only. The three sub-images show the cosine similarity matrices of the cross entropy trained (a), cosine embedding trained (b), and cosine embedding trained on the final ReLU layer (c) models. Since the matrices are normally symmetric, the bottom left half is set to the minimum value and crossed out.}
\label{fig:9}
\end{figure}

\subsection{Binary and Five-Way Classification Using Our Algorithm}
Since mere vector direction comparison seemed to work effectively with binary classification based on a cosine similarity matrix optimized towards one class (i.e. Figure 4), we expand the algorithm to include multi-class classification to use with the output vectors of AlexNet when trained on embeddings. Given Figure 9c, which shows the cosine similarity matrix obtained from the output vectors from this instance of AlexNet on the test set, the results should be notably better than the binary results for the engineered features (Figure 6). Now vectors in the same class, on average, are relatively aligned compared to vectors in different classes which have a comparably large angle between them.

From our previous cosine similarity matrices obtained with engineered feature vectors and the traditionally trained cross entropy and cosine similarity models, there had yet to be a cosine similarity matrix indicating a strong distinction between kink and column images. Figure 9c displays a strong distinction of kink and column images. This is important since there is a key physical distinction between kinks and columns which depends on the ratio of the current in the plasma divided by the flux of the background poloidal magnetic field through the an inner plasma gun electrode \citep[see eq. 2.17]{bellan}. As an example of how the algorithm should perform better given a more well behaved set of vectors, we perform kink versus non-kink classification using the binary algorithm from section 2.3 with $\gamma$ set to $100$ for optimal results and 1000 randomly sampled training pairs. We display the results via confusion matrix in the left of Figure 10. Given a kink image, the algorithm notably classifies it as a kink with 92\% accuracy. For the remaining classes, the algorithm correctly classifies images in our test set as non-kink with 96\% average accuracy. It should also be noted that in distinguishing between stable columns and columns exhibiting kink instability, the algorithm is 93\% accurate. With ten training samples instead of 1000, the accuracies are robust at 92\% and 97\% for kink and non-kinks, respectively.

Since a binary classification of vectors with well behaved cosine similarity achieved notable results, we modify the binary classification algorithm to perform multi-class classification based on an image being classified to the class it is most similar with. The modified algorithm is shown in the Algorithm 2 block. To change the algorithm for multi-class classification, we compute the average cosine similarity of one image with 50 samples from the training sets of each class to get five corresponding test means, $\mu_i$. We classify an image in the test set by assigning it to the class with the highest $\mu_i$. The results for full five-way classification of the test data (total of 285 images classified) are displayed in the right half of Figure 10 through a confusion matrix. The average accuracy for full, five-way classification of our test set is 92\%, roughly 10\% greater than our traditional transfer learning results of 82\% and 84\%. The results are identical when we compute the test means $\mu_i$ with 10 training samples from each class instead of 50.

\begin{algorithm}
\caption{: \textbf{Multi-Class Classification With Cosine Similarity} \\ Given: a number of classes $I$, a set of classes $C = \{C_i\}_{i \in I}$, a training set for each class $U_i$, a test set $P$, and a vector-valued function $f(x)$}\label{alg:2}
\begin{algorithmic}

\State \:\:\:\:\:\:\:\:1.\:\:\:\:\textbf{for} $p$ in $P$:
\State \:\:\:\:\:\:\:\:2.\:\:\:\:\:\:\:\:\:\:\:\:\textbf{for} $i$ in $I$:
\State \:\:\:\:\:\:\:\:3.\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\textbf{for} $u$ in $U_i$:
\State
\State \:\:\:\:\:\:\:\:4.\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:$x_u \xleftarrow{} \frac{f(p)\cdot f(u)}{\lVert f(p) \rVert \lVert f(u) \rVert}$
\State
\State \:\:\:\:\:\:\:\:5.\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:$\mu_i = \frac{1}{\lVert U_i \rVert}\sum\limits_{u=1}^{\lVert U_i \rVert} x_u$
\State
\State \:\:\:\:\:\:\:\:6.\:\:\:\:\:\:\:\:\:\:\:\: $p \in C_i$; where $i = \underset{i}{\mathrm{argmax}}\:(\mu_i)$
\end{algorithmic}
\end{algorithm}

\begin{figure}
  \centering
  \includegraphics[scale = 1.0]{fig10.png
  \caption{Confusion matrices of the test results using both the binary classification algorithm (left) and multi-class modification to the binary classification algorithm (right) on feature vectors obtained from AlexNet when trained on vector embeddings instead of the final layer. For the kink versus non-kink binary classification matrix the fractional accuracy for each class is displayed in the first row and the fractional inaccuracy for each class is presented in the second row. For the five-way classification results the diagonal values show the fractional accuracy of each class, while the off diagonal values show the fractional inaccuracies, displaying the fraction of each class in the x axis that was classified incorrectly as a class in the y-axis. The total amount of images classified for each case is 285.}
\label{fig:10}
\end{figure}

\subsection{Applying Previous Methods to the Large Data Set}

Our large data set from section 2.1 contains about 45000 unlabeled images. The instance of AlexNet trained using the cosine embedding loss on the final ReLU layer from section 2.4 is used to obtain feature vectors for each image in the large data set, and we classify the data set using the multi-class algorithm from section 2.5. We note the split of the data as 8\% haze, 15\% spider, 10\% column,  24\% kink, and 43\% sphere after classification is performed.

With feature vectors for each image, we represent a significant portion of the data in two dimensions using T-distributed Stochastic Neighbor Embedding (tSNE) \citep{tsne}, which is shown in Figure 11. We note that the classes are joined in a common center point which could be representative of the fact each class is almost orthogonal to each other class in the feature space, and by using the vectors alone we do not get a representation based on the actual purpose of the network. Thus, we also perform tSNE on the normalized feature vectors in order to give a better idea of the grouping based on direction alone, which is shown in Figure 11.
\begin{figure}
  \centering
  \includegraphics[scale = 1.0]{fig11.png
  \caption{T-distributed Stochastic Neighbor Embedding (tSNE) for 2000 sampled images from each class (left). The tSNE parameters of perplexity and verobosity are set to 350 and 2, respectively. We note that when the 4096 dimensional vectors are reduced to two dimensions the data appears to be generally separated appropriately. We conduct the same tSNE with normalized feature vectors, where instead of five lines of classes, we see five clusters of classes. The classes are labeled as 0, 1, 2, 3, 4 corresponding to the labeling convention from the previous sections in this work.}
\label{fig:11}
\end{figure}

In order to investigate the distribution of data from each class, we plot the distribution of vectors about the average vector for a given number of samples in each class. An example of this can be seen in left half of Figure 12, for the kink class. We note that the distributions all heavily skew towards an angle of zero degrees, meaning the vectors for each class mostly populate a small angular region, with few outliers relatively misaligned with the class. In order to get an idea for the quality of accuracy for a class, despite not having ground truth, we display the closest images to the average vector, randomly sampled images from zero degrees to ten degrees, and randomly sampled images from forty five to fifty five degrees in the right half of Figure 12 for the kink class. Given the network was trained to align vectors in the same class, these qualitative results suggest the images with vectors further away from the average vector should become more frequently classified incorrectly or be images that fall outside of the five predefined class.
\begin{figure}
  \centering
  \includegraphics[scale = 1.0]{fig12.png
  \caption{Distribution of 2000 randomly sampled vectors in the kink class about the average vector of the samples (left). The closest five images to the average vector, images randomly sampled from zero to ten degrees, and images randomly sampled from twenty five to thirty five degrees are displayed (right). 75\% of images in the class are within the zero to ten degree range, while the remainder of the data lies outside of this range.}
\label{fig:12}
\end{figure}

\section{Discussion}
This study is motivated by the goal of classifying a data set of plasma images from the Caltech Spheromak Experiment through feature extraction, transfer learning, and a vector direction comparison algorithm. We demonstrate the use of cosine similarity to select relevant features for binary classification with a well behaved cosine similarity matrix, which displays a representative average cosine similarity of vectors between each pair of classes. While the engineered features did not yield cosine similarity matrices indicating distinction among all five classes in our data set, simple vector direction comparison with no gradient descent based training was able to perform sufficiently well in the spider versus non-spider binary classification case. The idea of an extremely well behaved cosine similarity matrix inspired the use of AlexNet to train and yield one itself. We then use the embedding-trained model to conduct binary classification more accurately than the use of engineered features by 10\% from 86\% to 96\% simply because the feature vectors from this version of AlexNet were more well behaved. Based on the cosine similarity matrix in Figure 9c, when AlexNet is trained by the using the cosine embedding loss on the ReLU layer embeddings rather than the output layer, it seemingly groups image vectors in the same class together in the output space, so much so that representative average cosine similarity between test images in the same class is on average 0.84, corresponding to a relatively thin grouping for each class. Meanwhile the classes are separated enough for no class to share a similarity greater than 0.15 with another class. Thus, the image vectors from the final ReLU layer are divided into classes where each class is greatly separated by an angle relatively large compared to the angle between vectors in the class. This was exploited to expand our binary classification algorithm to a multi-class algorithm, and the five-way classification results are better than our traditional transfer learning results by about 10 percent from 82\% and 84\% to 92\%. When applying this uniquely trained version of AlexNet to the unlabeled data set, we can still see a propensity of the network to organize classes in thin groupings. In the ideal case with ground truth labels for a larger data set, the certainty an object belongs in the class as a function of angle away from the average vector can be probed quantitatively. This can potentially be used in the future to probe for sub-classes and outlier groupings corresponding to dense areas outside of the peak of a distribution. Preliminary developments of this idea are summarized in Fig. 11.

While the method has been tested by our data sets, a possible new direction could be changing activation functions at the end of AlexNet to allow the outputs to access a larger domain, the use with more CNNs and different truncations of them, and the extension of the method for testing with more than five classes (and larger data sets in general). In our case using AlexNet and the final ReLU layer was sufficient for five classes.

\section{Summary}
The work performed in this study focuses on the classification of images from the Caltech Spheromak Experiment using a feature vector direction comparison algorithm and a modification of traditional transfer learning with AlexNet. With two engineered features, we demonstrate the vector direction comparison algorithm to perform binary classification with 86\% accuracy. We improve the accuracy of the algorithm to 96\% by using vectors from an instance of Alexnet trained to align vector embeddings in the same class while separating vector embeddings among different classes. Full five-way classification of the data set is performed with the output vectors from the embedding-trained AlexNet and a modification to the binary algorithm, which achieves a test accuracy of 92\%. This was larger than our test accuracies with transfer learning using the cross entropy loss (84\%) and the cosine embedding loss (82\%) on the final layer of the network for traditional five-way classification.
\section*{Funding}
M. F. was supported in part by a U.S. DoE Science Undergraduate Laboratory Internships (SULI) award. This work was also supported in part by the U.S. Department of Energy through the Los Alamos National Laboratory (ICF program). Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).
\section*{Conflicts of Interest}
The authors declare no conflict of interest.


\section{Introduction}
A front is a transitional region between two states.
They can be infinitesimal in thickness or may take
the form of a wide gradient.
They may also be either static or dynamic,
 depending upon the nature
 of the system
 and states being modelled.  The study of the dynamics of fronts
  is important to many areas
  of sciences, in particular physics
  and applied mathematics,
  including the study of the effects of ionising radiation, of chemical catalysis by radiation, of the spread of fires, and of pressure shock waves.

A number of algorithms
   have been developed
   to model the propagation of fronts
   in various scenarios, many of which rely upon finding solutions to the Eikonal equation of motion.
Examples of Eikonal based models
  are the level-set method
  and its variants
  \cite{Sethian1996,GIBOU201882,doi:10.1137/0524066,Osher2003},
  which model the propagation of fronts
  from a set of initial conditions.
A simple yet useful class of algorithms
  in the restricted level-set methods
  are known as the fast-marching algorithms
  \cite{doi:10.1073/pnas.93.4.1591,10.2307/2653069},
  in which
  a monotonically evolving front
  is defined only by its local conditions, again through treatment of the Eikonal equation.
This simplification
  greatly reduces the computational cost.
A variant of the fast marching algorithms,
  which attempts to retain
  their computational advantages
  over other level set methods
  but allow for non-monotonically evolving fronts \cite{doi:10.1137/15M1017302}
  was also developed
  to give the flexibility
  to account for non-uniform
  interaction around the fronts.

However, these models would not be idea for describing the propagation of fronts
  driven by linear radiation emitted from a spatially finite source.
Examples of scenarios such models
  could describe is the propagation of ionisation fronts around radiative astrophysical objects such as stars
  \citep[see e.g.][]{Goldsworthy1961RSPTA}
  or photochemical reactions
  in a medium \citep[e.g.][]{Ducharme1990PhysRevB}
  irradiated by a light source.

In this paper,
  we present a new model
  of front propagation
  which was specifically
  developed to describe hydrogen ionisation
  around radiative astrophysical objects,
  but which can be applied to any scenario
  in which a source drives a front
  radially outwards.
This model divides the space
  around the object
  into a polar grid of cells,
  each of which contains a number of particles and,
  over a series of time steps,
  models the progression
  of the front as it transmits through these cells.
The precision of the model
  ultimately depends upon the sizes of the cells
  and time steps which,
  when the model is applied numerically,
  will be limited only by the computational resources available.
A similar method
  was developed in \cite{Altay:2015lua}
  for the study of astrophysical ionisation
  in media of homogeneous densities,
  but their study differed
  in its aims and applications
  from ours.

In Sec~\ref{sec:model},
  we introduce the model and describe its workings in a general scenario.
In Sec.~\ref{sec:star},
  we apply the model to the case of
  a Population-III star during the cosmological Epoch of Reionisation (EoR)
  and compare the results we obtain
  to those from a common analytic approximation,
  showing that the results
  are not only consistent
  but provide us
  with richer information
  than those pre-existing methods.
Finally, in Sec.~\ref{sec:disc},
  we highlight the capabilities of this model
  and its potential for modelling front propagation
  in complex media
  that researchers encounter
  in a variety of disciplines
  both within and beyond physics.

\section{The Radiative Front Model}
\label{sec:model}

Our model begins with a very generic assumption:
  radiation is being emitted from a point source
  and travelling radially away from that source.
As such, the model would be applicable
  to any system
  which features a source transmitting radiation into a surrounding medium which could be expected to illicit an alteration in that medium with minimal scattering.

The space around the source
  is divided into cells arranged into a polar grid,
  with the source at the origin.
Taking the example of such a grid in 3D space, the volume of each cell is given by
\begin{equation}
    V_{i,\theta,\phi}=\frac{4}{3}\frac{\pi}{N_{\theta}N_{\phi}}
    \left[\;\!
    (r_{i+1,\theta,\phi})^{3}-(r_{i,\theta,\phi})^{3}
    \;\! \right]
    \ ,
    \label{eq:cellvolume}
\end{equation}
   where $N_{\theta}$ is the number of cells
      defined in the $\theta$ dimension,
    $N_{\phi}$ is the number of cells
    defined in the $\phi$ dimension,
   and $r_{i,\theta,\phi}$ is the inner radius
      of the cell in question.
The medium under consideration is divided
  into particles
  distributed between the cells
    in accordance with the nature of the system being modelled.
The number of particles in a given cell
  with coordinates $\left(t,i,\theta,\phi\right)$
  that are in the initial state
  will be labelled $n_{\mathrm{I};t,i,\theta,\phi}$,
  and the number in the affected state
  will be labelled $n_{\mathrm{II};t,i,\theta,\phi}$.

Once the initial grid has been established,
  the radiating object is treated as being activated.
Assuming the speed of sound of its radiation
 is a known function, which we label as $v$,
 the maximum distance
 that radiation could have travelled radially
 from the source is computed at each time step
 as $v(t-1)\Delta t$, where $\Delta t$ is the length of each time step, such that $(t-1)\Delta t$ is the time that has elapsed
 by the beginning of that time step.
Every cell that is within this distance
 is then treated as being
 in contact with the radiation.

The medium in question has
  an absorbativity parameter associated with it.
In the case of a physical or chemical reaction
  this will be a function,
  $f$, of the particle cross-section
  of the substance and the particle number in a given cell.
We can therefore calculate the probability
   of each unit of radiation
   affecting a given unit of medium in a given cell as
\begin{equation}
    \Pi_{t,i,\theta,\phi}=1-e^{f_{t,i,\theta,\phi}} \ ,
\end{equation}
   such that the number of affected medium particles
   per time step will be given by
\begin{equation}
    \Gamma_{t,i,\theta,\phi}
    =J_{t,i,\theta,\phi}\;\! \Pi_{t,i,\theta,\phi}  \ ,
\end{equation}
    where $J_{t,i,\theta,\phi}$ is the flux of radiation entering that cell.

In order to calculate $J_{t,i,\theta,\phi}$,
   we will need to not only account for the angular distribution of radiation but also for the absorption
   in cells radially preceding the cell in question.
We define the radial absorption function as
\begin{equation}
    s_{t,i,\theta,\phi}=\sum_{j=1}^{i-1}f_{t,j,\theta,\phi}
\end{equation}
such that $s_{t,i,\theta,\phi}$ gives the fraction of radiation that reaches the $i$th radial cell in a given angular bin.  Thus, the flux of radiation entering the cell will be given by
\begin{equation}
J_{t,i,\theta,\phi}=J_{t,1,\theta,\phi}
 \ e^{s_{t,i,\theta,\phi}}
 \ .
\end{equation}

Having calculated the influx of radiation into a given cell at a given time step, we  recalculate the number of affected and unaffected particles as
\begin{align}
    n_{\mathrm{I};t,i,\theta,\phi}&=n_{\mathrm{I};,t-1,i,\theta,\phi}-\Gamma_{t,i,\theta,\phi} \ ,\\
    n_{\mathrm{II};t,i,\theta,\phi}&=n_{\mathrm{II};,t-1,i,\theta,\phi}+\Gamma_{t,i,\theta,\phi} \ ,
\end{align}
in a manner that the overall particle count in a given cell
  is preserved.
 Any other effects which might affect the particle counts,
   such as the stochastic processes
   that can reverse the radiation's effects,
   can be accounted for by adding additional terms into these equations.
Having recalculated the particle numbers in each cell
  at a given time step,
  the absorption fraction can either be calculated again as a function of that cell's parameters at that time step
  or more simply recalculated as
\begin{equation}
    f_{t+1,i,\theta,\phi}=f_{t,i,\theta,\phi} \
   \left(  \;\!
   \frac{n_{\mathrm{I};t,i,\theta,\phi}}{n_{\mathrm{I};0,i,\theta,\phi}} \;\!
   \right)  \ ,
\end{equation}
   accounting for the alteration in the number of reactive particles in the medium at a given time step by taking their ratio with the initial number,
   and the radial absorption function
   would then be recalculated in kind.

In effect, at every time step and every cell,
  we are calculating the amount of radiation reaching that cell
  from the source,
  taking into account absorption by preceding cells in the same angular bin, calculating the reaction rate in the cell due to that radiation, and preparing our map for a new set of calculations to take place at the next time step,
  taking into account the ensuing changes in the medium.
The mechanisms of this model is illustrated
  in the flow chart in Fig.~\ref{fig:flow}.

\begin{figure}[h]
    \centering
    \includegraphics{flowchart.pdf}
    \caption{A flow chart showing the basic underlying mechanism of the radiative front model.  $r$ is the radial distance to the inner edge of the cell in question, $t$ is the time step coordinate, $\Delta t$ is the real time duration of each time step, and $v$ is the speed of the radiation causing the front, such that the cell in question can be considered to be in contact with the radiation in any timestep for which $v(t-1)\Delta t>r$.  In the case of a background which allows for non-trivial transport of the medium through which the front is propagating, transport terms will be included in the cell parameter calculations to represent this.  In the event that the radiation is travelling at less than the speed of sound in that medium, an additional series of statements would need to be added within the loop over the map to find which cells are within range of medium transport and recalculate those cell's parameters even when they are out of range of the radiation.}
    \label{fig:flow}
\end{figure}

\section{Example Application - Ionisation in Astrophysical Media}
\label{sec:star}
\subsection{Background}
\label{sec:starbackground}

The earliest stars were composed almost entirely of hydrogen and helium,  and they are known as Population-III stars.
\cite{Glover2005SSRv,Wise2011nha}.
They were large and bright and emitted predominantly ultraviolet radiation.
They were also short lived, with lifespans generally
  not more than a few million years.
They are believed to have first come into existence
  during the transition
  from what are called the cosmological dark ages,
  during which the baryonic component of the Universe
  was primarily composed of neutral hydrogen (HI),
  to the EoR, during which that HI was predominantly converted into the ionised hydrogen (HII) that composes most of the baryonic matter of the Universe today
   \cite{Loeb:2000fc,Liu:2020krj}, beginning at a redshift\footnote{In cosmology, redshift is a measure of the expansion of the Universe; at a redshift of $z=15$, objects in the Universe will have been, on average, $1/15$th as far apart as they are now.} of approximately $z=15$ \cite{Koopmans:2015sua}.

An understanding of the ionisation dynamics
  around individual objects as Population-III stars
  is essential for high precision studies of the EoR,
  in particular, those being carried out
  by the 21cm surveys
  \citep{DiMatteo:2004jha,Pritchard2012RPPhy,refId0}
  though there is no precise model
  able to describe in detail
  the evolution, spatial extent and morphological structures
  of the fronts these sources create
  which determine the tomographic features
  to be observed in the 21cm spectral lines
  \cite{Mirocha:2017xxz}.

When ionising light travels through a medium containing neutral atoms,
  each photon will have a probability of ionising an atom of
\begin{equation}
    \Pi=1-e^{-\tau} \ ,
\end{equation}
where
\begin{equation}
    \tau=\sigma d
\end{equation}
is the optical depth of the path traversed for interaction cross-section $\sigma$ and the column density of the path $d=\rho_{\mathrm{HI}} l$
   for HI number density $\rho_{\mathrm{HI}}$ and distance $l$.
In the case of neutral hydrogen (HI)
  being converted to ionised hydrogen (HII),
  $\sigma=6.3\times 10^{-18}~\mathrm{cm^{-2}}$.

\begin{equation}
    \Gamma_{\mathrm{B}}=n_{\mathrm{e}}n_{\mathrm{p}}\alpha~,
    \label{eq:gammaB}
\end{equation}
for electron number density $n_{\mathrm{e}}$, proton number density $n_{\mathrm{p}}$, and B-recombination coefficient \cite{Osterbrock2006,1959MNRAS.119...90S,1964MNRAS.127..145P}
\begin{equation}
    \alpha\approx \frac{2.6\times 10^{-13}}{(T/10^{4}{\rm K})^{4/5}}
\end{equation}
with units of volume per unit time, where $T$ is the average temperature of the medium.

Due to B-recombinations, there will always be a theoretical limit to how large a photoionised bubble can be for a given radiation source.  That limit, given under the assumption that the system remains radiating and unchanging except for the interaction of the radiation and the medium for all of time, is given by the Str\"{o}mgren radius \cite{1939ApJ....89..526S}:
\begin{equation}
    r_{\mathrm{S}}=\left(\frac{3\phi}{4\pi n_{\mathrm{e}}n_{\mathrm{p}}\alpha}\right)^{\frac{1}{3}}~,
    \label{eq:RS}
\end{equation}
where $\phi$ is the total flux of photoionising particles being emitted.  This can be trivially rearranged to see that it is defined as the radius of a sphere around a radiative object at which the flux of radiation is equal to the number of B-recombinations; assuming that each photon can ionise only one particle, this will of course be the limiting radius of an ionisation bubble.  In the event that the ionising particles are able to ionise more than one atom each, this could be corrected by simply multiplying $\phi$ by the number of atoms each particle could ionise.

Of course, Eq.~\eqref{eq:RS} is not an entirely accurate model even for the idealised assumptions we are making; since the flux of particles will decrease with distance from the source, the number density of electrons and protons will also decrease, such that $r_{\mathrm{S}}$ should more realistically be phrased as an integral over spacetime with variable number densities.  Since Eq.~\eqref{eq:RS} would generally be calculated under the assumption that within the ionised region the number densities of ionised products equals the initial number density of neutral hydrogen and, if we were to proceed with such a more detailed analogue of the Str\"{o}mgren radius, we would expect that they would decrease with distance, we would find that such a calculation would yield a slightly higher result.  Therefore, we would expect the radius calculated with Eq.~\eqref{eq:RS} to actually be slightly smaller than the maximum radius reached by the ionisation front of a star.  Effectively, Eq.~\eqref{eq:RS} assumes a hard cutoff at which complete ionisation gives way to complete neutrality, while failing to take into account the existence of an intermediating front at all.

\subsection{Implementation}

We wish to model the physics described in Sec.~\ref{sec:starbackground} with the model we described in Sec.~\ref{sec:model}.  Specifically, we wish to model the evolution of the ionisation front around a single Population-III star in a neutral background medium during the EoR by treating the star as a source of radiation at the origin of a polar coordinate system and studying the effects of the radiation on the matter in the cell's that compose the spherical polar map.  The speed of the radiation is the speed of light such that cells are considered to be in contact with the source when $r_{i,\theta,\phi}\leq c(t-1)\Delta t$ where $c$ is the speed of light and we are assuming a negligible difference between the speed of light in a vacuum and the speed of sound of the light in the medium.

The absorption function of a given cell is given by its optical depth
\begin{equation}
    \tau_{t,i,\theta,\phi}=\sigma \;\!  \rho_{\mathrm{HI};t,i,\theta,\phi}l~,
\end{equation}
for cross-section $\sigma$ and cell depth $l=r_{i+1}-r_{i}$.  From the definition of $\tau$, we have that $e^{-\tau}$ is the fraction of radiation that transmits through the medium without scattering.

The radial absorption function is then given by a sum over the optical depths of the preceding cells in the same angular bin,
\begin{equation}
s_{t,i,\theta,\phi}=\sum_{j=1}^{i-1}\;\!
 \tau_{t,j,\theta,\phi} \ .
\end{equation}

In order to calculate the flux
  of radiation reaching each cell
  at each time step,
we must not only account for the division of the star's overall flux throughout the angular components of the map, but also for the absorption that takes place before they reach a given cell.  Thus, the radiation intensity at a given cell will be a function of the overall stellar luminosity\footnote{Luminosity is defined as being the overall radiation flux emitted by a light emitting source.} $L$, the number of angular components $N_{\theta}$ and $N_{\phi}$, and the radial optical depth $S_{t,i,\theta,\phi}$.

Thus, at a given cell which we label cell $i$, we have that the ionising radiation intensity is given by
\begin{equation}
    J_{t,i,\theta,\phi}=\frac{ Le^{-s_{t,i,\theta,\phi}}}{N_{\theta}N_{\phi}}\Delta t~.
\end{equation}
The ionisation rate within that cell is then given by
\begin{equation}
\Gamma_{t,i,\theta,\phi}=\frac{JN}{h\nu}(1-e^{-\tau_{t,i,\theta,\phi}}) \ ,
\label{eq:gammai}
\end{equation}
where $N$ is the number of atoms each photon is capable of ionising, given by $N=\lfloor h\nu/E_{\mathrm{i}}\rfloor$ where $E_{\mathrm{i}}$ is the energy needed to ionise a ground state atom.  In the event that multiple types of radiation are emitted, we can generalise Eq.~\eqref{eq:gammai} to become a sum over the different wavelengths emitted or, in the event of a continuous distribution of radiation wavelengths, this could be further generalised to
\begin{equation}
\Gamma_{t,i,\theta,\phi}
 =\int {\rm d}j
 \ \frac{J_{j}N_{j}}{h\nu_{j}}
 \;\! \left(
  1-e^{-\tau_{j,t,i,\theta,\phi}}
  \right) \ ,
\label{eq:gammai3}
\end{equation}
where $N_{j}=\lfloor h\nu_{j}/E_{\mathrm{i}}\rfloor$.  However, for this simplified study we content ourselves to treat the star's radiation as being purely composed of $13.6\mathrm{eV}$ ultraviolet light; given the heavily ultraviolet spectrum of a Population-III star, this is not a particularly inaccurate treatment and will allow accurate results to a level of precision that is satisfactory for the present study.

Given that HI has only one electron, we treat $n_{\mathrm{e}}=n_{\mathrm{p}}=n_{\mathrm{HII}}$, such that Eq.~\eqref{eq:gammaB} can be implemented as
\begin{equation}
    \Gamma_{\mathrm{B};t,i,\theta,\phi}= {(n_{\mathrm{HII};t,i,\theta,\phi})}^{2}
    \;\! \alpha_{t,i,\theta,\phi}
    \;\! V_{i,\theta,\phi}
    \;\! \Delta t \ .
\end{equation}

As a first approximation,
  we treat the temperature of all cells as being constant,
  at $10^{4}\mathrm{K}$.
This is an appropriate temperature
  of such ionised media\footnote{The temperature will be heavily dependent on the chemical composition of the medium; since this paper concerns itself with demonstrating the validity of the proposed method rather than a precise study of astrophysical ionisation, we content ourselves to use a constant temperature of $T=10^{4}~\mathrm{K}$ in both our analytic calculations with Eq.~\eqref{eq:RS} and our implementation of the radiative front model.}.
Furthermore, we treat the astrophysical medium as being entirely composed of hydrogen, such that at all times we will have $n_{\mathrm{e}}=n_{\mathrm{p}}$ and
$N=\lfloor\;\! h\nu/13.6\;\!{\rm eV} \;\!\rfloor$.
In a detailed study relevant to 21cm surveys, we would want to use a more detailed analysis with variable temperatures and a medium containing helium, but for the purposes of demonstrating the radiative front model and its application to the evolution of photoionised fronts, these approximations are valid, particularly when one bears in mind that the Str\"{o}mgren radius given in Eq.~\eqref{eq:RS} assumes a purely hydrogen background and is itself a function of temperature, which we can also set to a constant $10^{4}~\mathrm{K}$.

Thus, at each time step we are altering each cell's HI and HII count by
\begin{align}
    n_{\mathrm{HI};t,i,\theta,\phi}&=n_{\mathrm{HI};t-1,i,\theta,\phi}-\Gamma_{t,i,\theta,\phi}+\Gamma_{\mathrm{B};t,i,\theta,\phi}~,\\
    n_{\mathrm{HII};t,i,\theta,\phi}&=n_{\mathrm{HII};t-1,i,\theta,\phi}+\Gamma_{t,i,\theta,\phi}-\Gamma_{\mathrm{B};t,i,\theta,\phi}~,
\end{align}
in a manner that naturally preserves the overall number of particles per cell, before recalculating the optical depths as
\begin{align}
    \tau_{t+1,i,\theta,\phi}&=\sigma \rho_{\mathrm{HI};t,i,\theta,\phi}l~,\\
    s_{t+1,i,\theta,\phi}&=\sum_{j=1}^{i-1}\tau_{t+1,j,\theta,\phi}~,
\end{align}
in preparation for the next time step.

\subsection{Application}

We demonstrate the model by taking the example of a $100~\mathrm{M_{\odot}}$ Population-III star in a homogeneous background composed entirely of HI with a number density of $\rho_{\mathrm{HI}}=5\times 10^{-6}\;\!(1+z)^{3}\mathrm{~cm^{-3}}$
 \cite{Crighton:2015pza}, where the bracketed term accounts for the reduction in density with the expansion of the Universe. Such a star would have a lifespan of approximately $2.7~\mathrm{Myr}$ and would emit
$2.4\times 10^{51} \mathrm{~eV~s}^{-1}$
of predominantly ultraviolet radiation \cite{Wise2011nha},
which for the sake of simplicity we treat as being purely $13.6~\mathrm{eV}$ light.  With these parameters, we generated homogeneous density maps with densities equal to the average densities at $z=15$, $z=10$, and $z=5$.

Before using the model to obtain new results, we would like to confirm that it is consistent with existing ones.
In order to do this, we allow the stars to exist for unrealistically long lifespans and compare the maximum radii of their ionised bubbles to the Str\"{o}mgren radii calculated from Eq.~\eqref{eq:RS}.
In Fig.~\ref{fig:zmax}, we plot the maximum extent of the front against the Str\"{o}mgren radii calculated from Eq.~\eqref{eq:RS} at the same three redshifts, where the maximum extent of the front is defined in the algorithm as being the inner radius of the cells at which $1-n_{\mathrm{HI};t,i,\theta,\phi}/n_{\mathrm{HI};0,i,\theta,\phi}>0$ and $1-n_{\mathrm{HI};t,i+1,\theta,\phi}/n_{\mathrm{HI};0,i+1,\theta,\phi}=0$. In these cases, we adjusted the run time of the model to account for the respective times it took for the front to reach an equilibrium at the different redshifts, changing the length of each time step accordingly, as well as altering the radius of each cell.  In the $z=15$ case, we used cell radii of $2.5~\mathrm{pc}$ and time steps of $450~\mathrm{years}$, in the $z=10$ case we used a cell radius of $5.2~\mathrm{pc}$ and time steps of $1000~\mathrm{years}$, and in the $z=5$ case we used cell radii of $18~\mathrm{pc}$ and time steps of $5000~\mathrm{years}$.  Due to the isotropy of the map, we used $N_{\theta}=N_{\phi}=1$ in order to maximise computational efficiency.

We can see that the curves end up slightly exceeding the Str\"{o}mgren radius, as expected.  We also see that in higher densities, ionisation fronts tend to reach this equilibrium position sooner; in the case of $z=5$, for example, it took approximately $500~\mathrm{Myr}$ for the front to stop expanding, while in the case of $z=15$ it too approximately $50~\mathrm{Myr}$.  This is in keeping with the notion that, in higher density, regions, the higher rate of B-recombinations will bring about a total radiation absorption in a smaller time frame as well as a smaller volume.

\begin{figure}[h]
    \centering
    \includegraphics[width=.49\linewidth]{z=15max.pdf}
    \includegraphics[width=.49\linewidth]{z=10max.pdf}
    \includegraphics[width=.49\linewidth]{z=5max.pdf}
    \caption{The maximum extent of the ionisation front, defined as being the radius of the cell at which $1-n_{\mathrm{HI},t,i,\theta,\phi}/n_{\mathrm{HI,0},t,i,\theta,\phi}>0$ and $1-n_{\mathrm{HI},t,i+1,\theta,\phi}/n_{\mathrm{HI,0},t,i+1,\theta,\phi}=0$, over periods of time far greater than the realistic lifespan of a Population-III star assuming a constant density as given by the average densities at $z=15$ (top left panel), $z=10$ (top right panel), and $z=5$ (bottom panel), showing that, as expected, the ionisation fronts slightly exceed the Str\"{o}mgren radii calculated from Eq.~\eqref{eq:RS}, plotted alongside the model's results as grey lines.}
    \label{fig:zmax}
\end{figure}

Having confirmed that the model is consistent with expected results, we may use it to generate new results by studying the composition of the photoionised region and the front that borders it in more detail.   We therefore set the stars to exist for their realistic lifespan of $2.7~\mathrm{Myr}$ and ran the model with time steps of $2000~\mathrm{years}$ and cell radii of, in the case of $z=15$, $17.14~\mathrm{pc}$, in the case of $z=10$, $8.57~\mathrm{pc}$, and in the case of $z=5$, $5~\mathrm{pc}$, once again setting $N_{\theta}=N_{\phi}=1$ for efficiency since we are assuming a homogeneous cosmological background.

Running the model thus, we obtain the results shown in Fig.~\ref{fig:contours}.  These figures show the evolution of the angularly symmetric front with time, defined through the ratio $1-n_{\mathrm{HI}}/n_{\mathrm{HI,0}}=n_{\mathrm{HII}}/n_{\mathrm{HII,0}}$, the fraction of the medium which has been ionised.  As can clearly be seen, photoionised bubbles in higher density media as are presented by the higher redshift samples evolve more slowly and with a thinner front.  Indeed, not only do the higher density samples present smaller bubbles and fronts, but the ratio of the width of the front to the width of the photoionised region is much lower, indicating that the radiation is able to penetrate partially ionised cells much less effectively due to their higher count of HI particles.  We also find that the front during the earlier phases of expansion is much thinner than this and only becomes thicker with a more apparent gradient as the ionised bubble becomes larger, which is what one might intuitively expect as the radiation becomes more dispersed and so begins to ionise the media it is in contact with at a lower rate.

\begin{figure}[h]
    \centering
    \includegraphics[width=.49\linewidth]{contourz=15.pdf}
    \includegraphics[width=.49\linewidth]{contourz=10.pdf}
    \includegraphics[width=.49\linewidth]{contourz=5.pdf}
    \caption{The evolution in time of the photoionisation front created by a Population-III star emitting $2.4\times 10^{51} \mathrm{~eV~s}^{-1}$ of purely $13.6~\mathrm{eV}$ ultraviolet light consistently from its creation to its destruction at an age of $2.7 \mathrm{Myr}$ at $z=15$ (top left (panel), $z=15$ (top right panel), and $z=5$ (bottom panel), under the approximation of a constant background density at each redshift.  The hydrogen to the outside of the front is almost entirely neutral, that inside it is almost entirely ionised with a small amount of constant recombination preventing full ionisation from ever occurring, and the region within the front is defined as being the region that has a reduced ionised to neutral hydrogen ratio $1-n_{\mathrm{p}}/n_{\mathrm{H}}$ of between $0$ and $1$.  As one would expect, the photoionised bubble evolves more quickly and to a larger size in eras with a lower matter density and the front at a given time is wider, as the lower optical depth of the medium allows a greater depth of space to be being ionised simultaneously.  Cells in which $n_{\mathrm{HI}}=n_{\mathrm{HI,0}}$ at a given time step are not plotted and therefore make up the region beyond the blue curve of the outer layer of the front.}
    \label{fig:contours}
\end{figure}

Taking another example, that of an inhomogeneous background density, once again with a supersonic ionisation front and the same stellar parameters, we define a density map in which the background is primarily defined as being the average density at $z=5$ with a region of overdensity at $\theta=3\pi/2$ which decreases exponentially around that centre until it vanishes at $\theta=0$ and $\theta=\pi$.  Allowing the star to live for the unrealistic timeframe examined for the results shown in Fig.~\ref{fig:zmax} and studying again the maximum extent of the front, we find that a slice in time and the homogeneous angular dimension gives a final ionisation front at $2.7~\mathrm{Myr}$ which varies with density, as shown in Fig.~\ref{fig:zmaxasymmetric}.  This serves to highlight the models ability to study supersonic radiative transport in inhomogeneous media, such as in a region of space containing nebulae and other sources of over- and underdensity.

\begin{figure}[h]
    \centering
    \includegraphics[width=.49\linewidth]{z=5maxasymmetric.pdf}
    \caption{A slice in $\phi$ of the maximum extent of the ionisation front, defined as being the radius of the cell at which $1-n_{\mathrm{HI};t,i,\theta,\phi}/n_{\mathrm{HI,0};t,i,\theta,\phi}>0$ and $1-n_{\mathrm{HI};t,i+1,\theta,\phi}/n_{\mathrm{HI,0};t,i+1,\theta,\phi}=0$, after a period of time of $500~\mathrm{Myr}$ with a density contrast defined at $r\geq 9~\mathrm{kPc}$ around $\theta=3\pi/2$ and exponentially decreasing around that angle before vanishing at $\theta=0$ and $\theta=\pi$.}
    \label{fig:zmaxasymmetric}
\end{figure}

\section{Discussion}
\label{sec:disc}
We have proposed a new method for modelling the evolution of fronts around radiative objects and have examined the example application of the ionisation bubble around a star in a neutral medium.  We have shown that this method agrees with existing methods for estimating basic parameters of the front while providing significantly more information than those methods, allowing us to precisely model its evolution in time as a function of all relevant physical parameters.  We propose that this method could find many uses in varied fields, from modelling the effects of ionising radiation as examined here, to modelling chemical photocatalysis, to studying more varied phenomenon such as the propagation of ant expeditions as they explore the territory around their nests.  We intend to continue the development of this model and use it for astrophysical applications, specifically using it to study ionising objects during the cosmological epoch of reionisation and how the information this model provides can be applied to upcoming studies of 21cm spectroscopy and we encourage members of other academic groups to apply this method with appropriate modifications and parameters to their own fields.

The current formulation of the Radiative Front Model relies upon a static medium through which the front is propagating; this works for a supersonic front, but would run into issues when describing subsonic fronts as they are much more likely to encounter alterations due to fluid transport caused by the change of state they herald.  Furthermore, the model assumes that no radiation is scattered to non-radiative angles.

The subsonic front issue could be corrected for by introducing a set of medium transport equations into the model which allow the particle number to change between cells in a manner that preserves the overall particle number over the set of cells.

The scattering of radiation could be accounted for in one of a number of ways; the simplest would be to incorporate a scalar field into the grid which affects each cell as a function of the conditions of its surrounding cells, effectively modelling the non-radiative scattering by taking an average of its expected effects and applying them directly to each cell.  Alternatively, the radiation could be modelled coming out of each cell in which non-radiative scattering is occurring directly before being mapped onto the initial polar coordinate system of cells; however, this method would be significantly more computationally expensive.

Another possible modification which could improve both the precision and computational efficiency of the model would involve dividing each cell into smaller subcells which could be selectively explored depending upon the circumstances; for example, regions well within or well outside the front could be treated as single units, while transitional front cells could be divided up to get a smoother description of the front as it evolves.  This could be performed by assigning regions that have been transitioned enough to be considered no longer a part of the front but a part of the altered region with a logical statement and then grouping all cells in contact with each other with that logical statement together as one larger cell, while taking all cells that are in contact with the radiation but are still considered to be a part of the front into smaller cells.  However, for the purposes of this paper which was intended to describe and demonstrate our new method, we restricted ourselves to the simpler example of using fixed cell sizes.

We leave it to a future paper to explore these potential extensions of the model in more detail and limit ourselves to the study of entirely radiative, supersonic fronts in this paper.