\section{INTRODUCTION}
\label{sec:introduction}

Robust and accurate state estimation for micro aerial vehicles (MAVs) is of crucial importance as these versatile robots are taking the place of people to fulfill complex and dangerous missions, such as industrial inspection, remote sensing, search and rescue.
Due to their flexibility and ability to minimize risks to humans, there is an increased demand for their ability to achieve perception and localization safely and reliably in more challenging environments, such as fire scenes, subterranean settings, as well as GPS-denied environments. To localize in such environments, visible-light camera sensors seem to be a suitable choice for MAVs because of low power consumption, light weight and affordable price. A number of promising visual odometry and visual-inertial odometry methods have been proposed in recent years. Nevertheless, due to the poorly illuminated conditions and airborne obscurants conditions such as dust, fog and smoke in real world, the data of visible-light cameras can significantly degrade, which makes them unreliable for motion estimation.
\begin{figure}[htbp]
	\centering
	\includegraphics[width=1.0\columnwidth]{experiment_result/first_page2.jpg}
	\caption{The perception abilities comparison of different sensors in the same smoke-filled environments. (a) and (b) show the real scene of light-smoke and dense-smoke from a visible-light camera. (c) shows the ThermalPoint tracking results from a thermal camera. Red point means that the feature has been tracked more than ten times. Otherwise, it will be marked as a blue point. (d) shows the trajectory comparison of TP-TIO (blue) and LOAM (red). The red and yellow stars represent the final position of LOAM and TP-TIO, respectively. It can be seen that our method is closer to the origin (xyz axes) while LOAM has a large pose estimation drift and point cloud misalignment. (e) shows the trajectory comparison in the z-direction. It can be seen that our method has a small drift in the z-direction.}
	\label{fig:first_page}
	\vspace{-3mm}
\end{figure}

Compared with visible-light cameras, thermal cameras are robust to illumination variations and obscurants like dust and fog. Thermal cameras, especially the thermal cameras operating in the long wave infrared spectrum (LWIR), have better spectral penetration properties in the presence of fog, dust, and smoke (shown in Fig.\ref{fig:first_page}). Therefore, thermal odometry is an ideal choice to achieve robust pose estimation in visually degraded environments.
In spite of this, the use cases of the thermal camera are largely limited to perception and localization. Since the thermal camera captures 16-bit temperature radiation images which have low contrast, it is difficult to use handcrafted feature methods to extract salient features. To solve this problem, most of the previous methods \cite{borges2016practical,mouats2014multispectral} re-scale 16-bit thermal images to 8-bit to enhance the contrast of thermal images and improve the performance of feature matching. However, re-scaling thermal image data will result in loss of information, amplifying the effect of noise and image saturation \cite{khattak2019keyframe}. Moreover, thermal images have fixed-pattern noise and requires a  Non-Uniformity Correction (NUC) process \cite{borges2016practical} to remove noise, which  periodically suspends the images for up to 500 milliseconds. This kind of image interruption can easily cause the failure of feature tracking and drift of state estimation.

In the last decade, due to the ability to extract high-level features automatically, deep neural architectures have been proved advantageous for feature extraction in various data types including RGB images and point clouds.
However, the shortcomings of deep learning based feature extraction networks are also obvious, such as high computation cost and poor real-time performance.

Motivated by the discussion above, we propose a novel deep thermal-inertial odometry (TP-TIO).
The main contributions of our work are as follows:

\begin{itemize}
	\item
	    We propose the first tightly coupled deep thermal-inertial odometry algorithm (TP-TIO), which utilizes both a lightweight learning-based feature detection network and an optimization-based visual-inertial framework.

	\item We propose the ThermalPoint, a CNN-based feature detection network specifically tailored for producing keypoints on thermal images, that provides notable improvements in anti-noise compared with  other state-of-art learning-based and handcrafted feature extraction methods.

	\item We propose an IMU-aided radiometric feature tracking method, which directly makes use of full radiometric data from a infrared camera to establish robust correspondences between consecutive frames.

	\item We conducted an extensive set of experiments to thoroughly evaluate our framework from thermal feature tracking to pose estimation in all kinds of environments including obscurant-filled underground mines and dense smoke-filled indoor environments.

\end{itemize}
The rest of the paper is organized as follows: Section \ref{sec:relatedwork}
reviews various feature extracton methods and thermal-inertial odometry methods. Section \ref{sec:algorithm} illustrates the details of the proposed
method. Quantitative and qualitative evaluations from ThermalPoint to state estimation are presented in Section \ref{sec:experiment}. Section \ref{sec:conclusion} concludes the paper.
\section {Related Work}
\label{sec:relatedwork}
In this section, we start by reviewing existing feature extractor methods, which are essential for vision-based and thermal-based motion estimation, then we review various thermal odometry or thermal SLAM methods.
\subsection{ Feature Extractor on Thermal Images}
In decades, various handcrafted feature point detection and descriptor algorithms have been proposed and successfully utilized in many visual odometry and SLAM framework. However, since the classical feature extractors are based on gradient information of the images, they cannot work well when images have severe transforms (e.g. viewpoint, blurring, noise). To solve this problem, recent works\cite{superpoint, LIFT,gcnv2} leverage deep learning methods to generate more robust interest points and descriptors. Among these methods, LIFT \cite{LIFT} performs both keypoint detection and descriptor extraction based on patch methods and requires a supervision strategy. SuperPoint \cite{superpoint} proposes a homographic adaptation strategy and predicts the keypoints and descriptors from an unsupervised approach. GCNv2 \cite{gcnv2} adopts a binary descriptor with good computational efficiency. However, all of these methods are designed for standard visible band images and they are difficult to achieve good performance on thermal imagery.

\subsection{Thermal Odometry}
Robust and accurate state estimation from a thermal camera remains a challenging problem. In recent years, there are a number of efforts that have been made for the thermal odometry.
To improve the robustness of thermal odometry, many works fuse thermal images with other modalities like RGB images\cite{mouats2014multispectral,poujol2016visible,dai2019multi} or laser scanner\cite{shin2019sparse}.
Since the above methods use traditional handcrafted feature detection methods that mainly rely on the gradient of images, it is difficult for these methods to overcome the effect of large noise. To solve this problem, we propose the "ThermalPoint", which adds a noise reduction module during the training process and effectively overcomes the noise of thermal images. Khattak et al.\cite{khattak2019keyframe} developed a keyframe-based thermal inertial odometry and proved the effectiveness of the direct method. However, since the direct methods heavily rely on image quality and accurate initial pose from the back-end module, it is particularly challenging for these methods to provide robust motion estimation when the environments have a scarcity of thermal gradients. In contrast, since our feature tracking method doesn't rely on the feedback of the back-end module, it can handle temporary feature tracking failure in environments with no salient heat sources.
Muhamad\cite{saputra2020deeptio} proposed an end-to-end thermal inertial odometry. However, this method can only perform well with a slow frame rate. To achieve motion estimation at a high frequency, our method uses ThermalPoint in a shallow architecture and radiometric feature tracking method to establish reliable correspondences instead of matching descriptors. This greatly improves real-time performance and reduces computation.

In summary, it is relatively difficult to find a real-time and robust deep feature-based thermal-inertial odometry under severely visually degraded conditions in the existing literature. This paper aims to solve this problem, providing a reliable sparse feature-based state estimation method in the thermal image domain.

\section{Thermal-Inertial Odometry} \label{sec:algorithm}

\begin{figure*}
\centering
\includegraphics[width=2.0\columnwidth]{experiment_result/framework/thermal_odometry.jpg}
\caption{Overview of deep thermal-inertial odometry. The modules marked as yellow are our contributions to this work. To improve real-time performance, the modules with green dashed squares will only participate in the training process and be excluded in a real application.}
\vspace{-3mm}
\label{fig:algorithm}
\end{figure*}

In this section, we will introduce the pipeline of robust Thermal-Inertial Odometry with deep "ThermalPoint" (TP-TIO) as shown in Fig.\ref{fig:algorithm}, which allows robust motion estimation in severely visually degraded environments. It is composed of three sequential modules, namely ThermalPoint network, IMU-aided radiometric feature tracking and thermal-inertial bundle adjustment.
\subsection{ThermalPoint Network}

\subsubsection{Network Architecture}

To detect features robustly on noisy thermal images, we designed a fully convolutional neural network architecture called ThermalPoint shown in Fig.\ref{fig:algorithm}. Inspired by the famous SuperPoint \cite{superpoint} network, this network has a shared encoder which abstracts high-level features and reduces the dimensionality from the input image. After that, the architecture is divided into two task decoders, feature detection decoder and descriptor decoder, respectively. To improve real-time performance,
the descriptor decoder is only used to improve the repeatability\cite{superpoint} of feature detection during the training process and is removed in the inference process. Instead, we use KLT-based feature tracking method to establish feature correspondences. More details can be found in Section \ref{sec:tracking}.

\subsubsection{Shared Encoder}
Considering the low latency and computation cost for feature detection on embedded systems, our shared encoder is built from an efficient network, MobileNet\cite{howard2017mobilenets}, instead of the  VGG-style shared encoder used by SuperPoint. The specific structure of our shared encoder can be found in Fig.\ref{fig:algorithm}.

\subsubsection{Feature Detection Decoder And Descriptor Decoder}

Inspired by SuperPoint \cite{superpoint}, the feature detection decoder computes a feature map $X\in \mathbb{R}^{H_c\times W_c \times 64}$.
 After that, it will output the final confidence heatmap $C$ sized $H \times W$ by performing a channel-wise softmax and DepthToSpace. Each pixel in $X$ corresponds to a probability, measuring a likelihood to be the feature point in the $8\times 8$ cell. For the descriptor part, ThermalPoint is similar to the SuperPoint. Here are some differences we have improved between ThermalPoint and SuperPoint \cite{superpoint}.

\begin{itemize}
	\item
    To improve the repeatability\cite{superpoint} of feature points, we replace multiple cross-entropy loss used by SuperPoint with an independent logistic loss (Eq.(\ref{eq:detector loss})). The reason is that multiple cross-entropy loss only allows one feature extracted in the cell, while independent logistic loss encodes the probability for every pixel.
	\item To save computation time during training, our decoder performs lower computational bilinear interpolation $Y$ to get dense descriptor map $D$, instead of bi-cubic interpolation of the SuperPoint.
\end{itemize}

\subsubsection{Loss Function}
Similar to the SuperPoint's training process, we optimize both the feature point detector loss and the descriptor loss.
For the descriptor loss, we directly employ the SuperPoint's loss. The feature point detector loss function $L_{dec}(I)$ is different, where we use a binary logistic loss over each pixel of $C$. Given the feature point label $L(I) \in \mathbb{R}^{H \times W}$for each image $I$, the loss can be defined as:
\begin{equation}
\begin{aligned}
\label{eq:detector loss}
L_{dec}(I) = \frac{1}{HW} \sum_{i=1,j=1}^{H,W} -log(C(i,j))L(i,j) - \\
log(1 - C(i,j))(1-L(i,j))
\end{aligned}
\end{equation}

\subsection{Training Details}
\subsubsection{Noise Reduction} \label{Data_augmentation}

 Applying various noise models to training datasets is a common technique to improve the model's anti-noise performance. Compared with the noise presented in visual images, thermal images have more complex noises, especially Fixed Pattern Noise (FPN) which is also called non-uniform noise. To overcome the effects of FPN on feature extraction, we need to add FPN noise to the training datasets (Fig.\ref{fig:fpn}). Given the non-uniform image $\tilde{I}(i,j)$ and the noise-free image $I(i, j)$, the relation of them can be simply described as the below equation
\begin{equation}
\label{eq:nonuniform}
\tilde{I}(i,j) = I(i,j) + n(i,j)
\end{equation}
where $n(i,j)$ is the non-uniform noise.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=.4\textwidth ]{experiment_result/data_augmentation/FPN.pdf}
	\caption{Samples of applying FPN on datasets. The first row shows original noise-free images from synthetic dataset, KITTI dataset and thermal dataset. The bottom row shows corresponding non-uniform images after adding FPN noise.}
	\label{fig:fpn}
\end{figure}

To obtain non-uniform noise, we employ two simple and effective methods: The first is facing a thermal camera toward a scene with the same material like a concrete wall. In terms of their heat capacity, the non-uniform noise will be recorded due to the scarcity of thermal gradients. The second is to simulate column FPN. In this paper, we use both the above two methods to create a non-uniform noise dataset which contains 1000 images. During the training process, we randomly select an image from the non-uniform noise dataset and add it to a noise-free image based on Eq.(\ref{eq:nonuniform}) (See Fig.\ref{fig:fpn}). Meanwhile, we also use standard data augmentation tricks including contrast, hue, brightness to improve the network’s robustness against photometric changes.
\subsubsection{Pseudo Feature Supervision} \label{multiple_feature_point_labels}
SuperPoint \cite{superpoint} shows great benefits of using anchors to guide its training process. Although anchors provide a stabler training process, they also prevent the model from identifying new keypoints when no anchor is in the proximity. To solve this problem, we use SIFT \cite{SIFT} to generate some extra pseudo keypoint labels and supervise the training process. By combining the labels from handcrafted features and learned CNN filters, the repeatability of ThermalPoint is improved, as shown in Table\ref{tab: feature detectors repeatability}. The specific detail can be found in Fig.{\ref{fig:label}}.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=.48\textwidth ]{experiment_result/data_augmentation/label.jpg}
	\caption{Pseudo feature supervision training on RGB (left) and Thermal dataset (right), respectively. The green and red points represent keypoints detected by our base feature detector and SIFT \cite{SIFT} method, respectively. Our feature detector mainly detects some "conner" points, while the SIFT method can detect some "salient gradient" points. Both  keypoints will be used as labels during the training process.}
	\label{fig:label}
\end{figure}

	\subsection{IMU-Aided Radiometric Feature Tracking}
	\label{sec:tracking}

After detecting a sparse set of keypoints through our thermal point network, we need to track the motion and build the correspondence of these points between frames. Since most of deep learning-based feature detection networks \cite{LIFT, superpoint} create dense descriptor map to find correspondences of features, it is difficult to meet our real-time requirements. Inspired by the work\cite{khattak2019keyframe}, instead of using the descriptor method, we use full radiometric based KLT method and IMU measurement to achieve robust feature tracking. To overcome the sudden photometric changes, we combine the inverse-compositional approach \cite{inversecomposition} with a patch dissimilarity norm, making our feature tracking invariant to intensity scaling.

    Given a set of patches in a reference image, we minimize the radiometric error between corresponding patches in two consecutive frames as estimating the transform $\mathbf{T} \in \mathrm{SE}(2)$. For a single residual in a patch, it can be defined as

	\begin{equation}
	r_{i}(\boldsymbol{\xi})=\frac{I_{t+1}\left(\mathbf{T} \mathbf{x}_{i}\right)}{\overline{I_{t+1}}}-\frac{I_{t}\left(\mathbf{x}_{i}\right)}{\overline{I_{t}}} \quad \forall  \mathbf{x}_{i} \in \Omega
	\end{equation}

	Where $I_{t}(\mathbf{x})$ is the thermal value of the 16-bit image $t$ at pixel location $x$. The mean thermal value if the patch in image t is $\overline{I_{t}}$. $ \Omega$ defines the set of image coordinates in the patch.

	To improve the robustness for large displacements and real-time performance, we perform our patch alignment over an image scale-space pyramid from coarsest level to finest level, where the alignment results from the upper level provide a prior for the lower level. Furthermore, we use the IMU preintegration \cite{VINS} to provide transformation prior for the alignment process at the coarsest level. These two steps greatly improve the solution convergence speed and improve the robustness of the tracking.

	\subsection{Thermal-Inertial Bundle Adjustment}

	To estimate the motion of the sensor, an optimization is performed that jointly minimizes visual reprojection error terms of the tracked features and IMU error terms based on IMU preintegration measurements within a sliding window.

	\subsubsection{Reprojection error} The first residual term which we use for pose estimation is the reprojection error based on thermal images. Thermal reprojection error can be described as
	\begin{equation}
	\mathbf{e}_{\text {reproj }}=p^{{c_{j}}}-K\left(T_{{b}}^{c} T_{{w}}^{b} l^{{w_{i}}}\right)
	\end{equation}
	Where $p^{{c_{j}}}$ is the point in the target frame $j$ at camera image coordinates $c$. $l^{{w}}$ is the corresponding 3D landmark location in  the reference frame $i$ at world coordinates $w$. $T_{{w}}^{b} $ is the transformation from world frame $w$ to IMU body frame $b$. $T_{{b}}^{c}$ is the transformation from camera frame $c$ to IMU body frame $b$.

	\subsubsection{IMU error} The second residual term is derived from IMU data based on IMU preintegration measurements. To achieve high frequency and robust pose estimation, we preintegrate several consecutive IMU measurements and add the IMU factor between frame $i$ and frame $j$. Similar to \cite{VINS}, the corresponding preintegration between IMU body frame $b_{k}$ and $b_{k+1}$ is $\Delta \mathbf{r}=(\hat{\boldsymbol{\alpha}}_{b_{k+1}}^{b_{k}}, \hat{\boldsymbol{\beta}}_{b_{k+1}}^{b_{k}}, \hat{{\gamma}}_{b_{k+1}}^{b_{k}})$. Therefore, the error term for inertial measurement within two consecutive frames can be defined as:

	\begin{equation}
	\mathbf{e}_{\text {imu }}=\left[\begin{array}{c}
	{\delta \boldsymbol{\alpha}_{b_{k+1}}^{b_{k}}} \\
	{\delta \boldsymbol{\beta}_{b_{k+1}}^{b_{k}}} \\
	{\delta \boldsymbol{\theta}_{b_{k+1}}^{b_{k}}} \\
	{\delta \mathbf{b}_{a}} \\
	{\delta \mathbf{b}_{g}}
	\end{array}\right]
	\end{equation}
	$=\left[\begin{array}{c}
	{\mathbf{R}_{w}^{b_{k}}\left(\mathbf{p}_{b_{k+1}}^{w}-\mathbf{p}_{b_{k}}^{w}+\frac{1}{2} \mathbf{g}^{w} \Delta t_{k}^{2}-\mathbf{v}_{b_{k}}^{w} \Delta t_{k}\right)-\hat{\boldsymbol{\alpha}}_{b_{k+1}}^{b_{k}}} \\
	{\mathbf{R}_{w}^{b_{k}}\left(\mathbf{v}_{b_{k+1}}^{w}+\mathbf{g}^{w} \Delta t_{k}-\mathbf{v}_{b_{k}}^{w}\right)-\hat{\boldsymbol{\beta}}_{b_{k+1}}^{b_{k}}} \\
	{2\left[\mathbf{q}_{b_{k}}^{w_{k}^{-1}} \otimes \mathbf{q}_{b_{k+1}}^{w_{w}} \otimes\left(\hat{\gamma}_{b_{k+1}}^{b_{k}}\right)^{-1}\right]_{x y z}} \\
	{\mathbf{b}_{a b_{k+1}}-\mathbf{b}_{a b_{k}}} \\
	{\mathbf{b}_{w b_{k+1}}-\mathbf{b}_{w b_{k}}}
	\end{array}\right]
	\\
	\\
	$
	Where $(\cdot)^{w}$ and $(\cdot)^{b}$ are world frame and IMU body frame respectively. $[\cdot]_{x y z}$ extracts the vector part of quaternions $\mathbf{q} $. Both rotation matrices $\mathbf{R}$ and Hamilton quaternions $\mathbf{q} $ are used to represent rotation. $\mathbf{p}_{b_{k}}^{w}$, $\mathbf{v}_{b_{k}}^{w}$ and $\mathbf{q}_{b_{k}}^{w}$ are translation, velocity and rotation from the body frame to the world frame while taking the $k$th image. $\mathbf{g}^{w}$ is the gravity vector in the world frame.
	$\mathbf{b}_{a}$ and $\mathbf{b}_{w}$ are accelerometer bias and gyroscope bias respectively. For more specific details about the preintegration theory, we refer the reader to \cite{VINS}.

	\subsubsection{Joint Optimization}

Since we have obtained the reprojection and inertial residual equations, for each new frame, we minimize the cost function that consists of reprojection residuals, IMU residuals and a marginalization prior $\mathbf{e}_{m}$.

	\begin{equation}
	J=\sum_{ i \in \mathbf{obs}(i)}\mathbf{e}_{\text {reproj }}^{\mathbf{T}} W_{\text {reproj}}^{-1} \mathbf{e}_{\text {reproj }}^{\mathbf{}}+\sum_{(a, b) \in \mathbf{C}}\mathbf{e}_{\text {imu}}^{\mathbf{T}} W_{imu}^{-1} \mathbf{e}_{imu}^{\mathbf{}}
	+\mathbf{e}_{m}
	\end{equation}

	$obs(i)$ represents a set that contains the ThermalPoint $i$ that is tracked and observed by other frames. The set $\mathbf{C}$ contains pairs of frames $(a, b)$ connected by IMU factors. $W_{\text {reproj}}$ and $W_{\text {imu}} $ are co-variance matrices for visual and IMU measurements. Through optimizing this cost function, we can obtain the sensor pose.

	\begin{figure*}[hbt]
		\centering
		\includegraphics[width=1.95 \columnwidth, height=0.5 \columnwidth]
		{experiment_result/feature_illumi_test_jfr_mine/feature_illumi_jfr_mine.png}
		\caption{The anti-photometric comparison of different feature detectors with sudden photometric changes.}
		\label{fig:feaure_illumi}
	\end{figure*}

	%

	\section{Experiments}
	\label{sec:experiment}
	To thoroughly evaluate the proposed TP-TIO framework and validate the performance of ThermalPoint, a variety of experiments were performed including underground mine, indoor laboratory and dense smoke-filled environments.
	\begin{itemize}
		\item For the thermal feature detection and tracking, we mainly compare the robustness and real-time performance of ThermalPoint detection and tracking with other state-of-the-art methods including Harris \cite{harris}, Good feature \cite{goodfeature}, ORB \cite{ORB}, SIFT \cite{SIFT} and SuperPoint \cite{superpoint}.
		\item For the state estimation, we mainly compare the accuracy and robustness of TP-TIO with other state-of-the-art odometry methods such as VINS \cite{VINS}, OKVIS \cite{OKVIS}, R-VIO \cite{R-VIO} and DSO \cite{DSO}.
	\end{itemize}
	The software runs on Jetson AGX Xavier computer with a Linux system running Robot Operating System (ROS).
	Our experiment video is available at: \textcolor{blue}{\href{https://youtu.be/aa4whgmYTqY}{https://youtu.be/aa4whgmYTqY}}

	\subsection{Robustness Comparison of ThermalPoint Detection}

	To verify the robustness of ThermalPoint detection, we evaluate anti-photometric performance and anti-noise performance.
	For the anti-photometric test, we use a public underground mine dataset containing severe thermal photometric changes. For the anti-noise test, we use our indoor smoke dataset presenting the challenges of large noise and dense smoke. ThermalPoint is evaluated against other famous state-of-the-art feature detection methods such as Harris, GoodFeature, ORB, SIFT and SuperPoint. We select a suitable fixed threshold for each feature detector to ensure most features can be extracted on this dataset. For this test, we compute a maximum of 400 points for all systems at a 512x640 resolution in each image. The Non-Maximum Suppression(NMS) of all systems is 8.

	%

	\subsubsection{ Anti-photometric and Anti-noise Comparison} \label{Anti-photometric comparison}
	Fig.\ref{fig:feaure_illumi} shows the anti-photometric comparison of different feature detectors with sudden photometric changes. We can easily find that deep learning-based methods (ThermalPoint and SuperPoint) can robustly detect enough features even on the low contrast thermal images. However, the classical feature detectors like Harris, Good Feature, ORB and SIFT struggle in the presence of sudden photometric changes. Fig.\ref{fig:feaure_noise_test} shows the anti-noise comparison of different feature detectors in the presence of fixed-pattern noise appearing as vertical stripes. It can be seen that ThermalPoint extracts few features on the stripe while other methods detect many features on the stripe.

The reason behind this is that our deep-learning-based method can learn anti-photometric and anti-noise abilities through a training process. In contrast, most of the classic detectors are based on the gradient information of images; Therefore, it is difficult for classic detectors to overcome sudden photometric changes.

	\begin{figure}[h]
		\centering
		\includegraphics[width=.48\textwidth]{experiment_result/feature_noise_test_jfr/anti-noise-cmu.png}
		\caption{The anti-noise comparison of different feature detectors in the presence of fixed-pattern noise. The area marked by the red rectangle presents the effect of fixed-pattern noise, which appears as vertical stripes of different brightness in the thermal images.}
		\label{fig:feaure_noise_test}
	\end{figure}

	\subsubsection{Quantitative Robustness Analysis} \label{Quantitative Analysis of feature detection }
	To further quantitatively analyze the robustness of feature detectors under viewpoint changes and photometric changes on thermal images, we use the repeatability score \cite{superpoint} to evaluate the robustness for each method. If the repeatability score is high, the corresponding method has a good ability of anti-photometric and anti-viewpoint changes. To test the repeatability on thermal images, we apply random homography transformations on our outdoor thermal image sequences to obtain corresponding thermal image pairs. We compute a maximum of 500 points for all systems at a 240x320 resolution. The Non-Maximum Suppression of all systems is 8 and the correctness threshold \cite{superpoint} $\varepsilon$ is 3. Table \ref{tab: feature detectors repeatability} presents the repeatability score for each method with all kinds of changes (viewpoint and photometric). It can be seen that ThermalPoint is the most repeatable under photometric and viewpoint changes.

	\begin{table}[h]
		\begin{center}
			\caption{Repeatability results of different feature detectors on our simulated dataset }

			\label{tab: feature detectors repeatability}

			\setlength{\tabcolsep}{10pt}

			\begin{tabular}{llllllllllllp{1cm}}
				\toprule
				Approach  			& Repeatability            	     \\
				\midrule
				Harris              & 0.4313                         \\
				Good feature  		& 0.6253   	                     \\
				ORB  				& 0.5748                 	     \\
				SIFT  				& 0.4976   	                     \\
				SuperPoint  		& 0.7166   	                     \\
				Our method          &\textbf{0.7674}                \\
				\midrule
			\end{tabular}
			\vspace{-8mm}
		\end{center}
	\end{table}

	\subsection{Robustness and Run-time Comparison of ThermalPoint Tracking}
	\label{Analysis of feature tracking }

	To verify the robustness of ThermalPoint tracking, we test our method in an underground mine environment shown in Fig. \ref{fig:feaure_math_jfr}. During this test, we track a maximum of 160 points for all systems at a 640x512 resolution for every frame. The Non-Maximum Suppression (NMS) of all systems is 8. We find that ThermalPoint can track significantly more points than other algorithms on thermal images with low contrast. Also, we evaluated the real-time performance of various feature tracking algorithms (Fig. \ref{fig:feaure_track_time_statistic}). Our feature tracking algorithm takes significantly less time than SuperPoint.

	\begin{figure}[b]
		\centering
		\includegraphics[width=.48\textwidth]{experiment_result/feature_match/jfr_feature_match2.jpg}
		\caption{Samples of different feature tracking method in underground mine dataset.}
		\label{fig:feaure_math_jfr}
	\end{figure}

	\begin{figure}[h]
		\centering
		\includegraphics[width=.4\textwidth, height= 0.36 \columnwidth ]{experiment_result/feature_track_time_statistic/time_statistics.png}
		\caption{Average time for
detecting and tracking features between two frames}
		\label{fig:feaure_track_time_statistic}
	\end{figure}
	\subsection{Accuracy Comparison of State Estimation in Indoor Environments}

	\begin{figure}[htbp]
		\centering
		\includegraphics[width=.48\textwidth ]{experiment_result/odometry_comparison/accuracy_comparison.jpg}
		\caption{(a1) and (a4) show the accuracy comparison of estimated trajectories of Ours (blue line), VINS (green line), OKVIS (yellow line) and R-VIO (purple line) against ground-truth (red line) in our laboratory environment as shown in (a2). DSO trajectory was excluded from the plots due to multiple failures. (a3) shows the quality of thermal images and  the feature detection result of our method in that scene.}
		\label{fig:accuracy}
	\end{figure}
To provide a comprehensive accuracy evaluation, the proposed method was compared with four state-of-the-art visual odometry methods, namely VINS \cite{VINS}, OKVIS [25], R-VIO [23] and DSO [19] on thermal imagery. The selections of the above methods are made mainly for two reasons: 1) Most of the thermal odometry approaches \cite{mouats2014multispectral,poujol2016visible,dai2019multi} require both visual and thermal images to achieve good motion estimation. However, our method works purely on thermal imagery; 2) Since most of the above methods use classic feature detection techniques, by comparing these methods, the effect of feature detection on thermal odometry will be clearly presented.

	To achieve a more fair comparison for state-of-the-art methods, we choose our laboratory environment with more heated sources (Fig. \ref{fig:accuracy}(a2)) and obtain high contrast thermal images (Fig. \ref{fig:accuracy}(a3)), which will be easier for state-of-the-art methods to achieve feature detection and tracking on thermal images. We held our device and walked 85 m along the bridge (Fig. \ref{fig:accuracy}(a2)) and returned back to the same position including purely rotational motion. The ground-truth is  provided by LOAM \cite{loam} algorithm.

	Fig. \ref{fig:accuracy}(a1,a4) shows the accuracy comparison of estimated trajectories in indoor environments. We can easily find that the trajectory provided by our method (blue line) is the closest to the ground-truth (red line), which suggests that our method is more accurate than other methods.

	The Absolute Trajectory Error (ATE)  and RPE \cite{benchmark} are presented in Table \ref{tab:accuracy}. It can be noted that the overall ATE and RPE errors of our method are very small.

	%

	\begin{table}[h]
		\begin{center}
			\caption{Accuracy evaluations of OURS, VINS,OKVIS,R-VIO and DSO methods operating on thermal images in indoor environments}
			\setlength{\tabcolsep}{10pt}
			\newcommand{\tabincell}[2]{\begin{tabular}{@{}#1@{}}#2\end{tabular}}
			\begin{tabular}{llllllllllllp{1cm}}
				\toprule
				\multirow{2}{*}{} 						& \multicolumn{2}{c}{ATE (in m)} 	& \multicolumn{2}{c}{RPE (in m)}    	     \\ \midrule
				Approach  					  			& RMSE     	& MAX     	& RMSE          &  MAX    \\ \midrule
				OURS                				& \textbf{0.359}   	&\textbf{0.866}    & \textbf{0.007}  &  \textbf{0.153}     \\
				VINS  					& 2.406   	& 3.902   & {0.033}   &  0.165         \\
				OKVIS  					& 3.199   	& 5.517   & {0.011}   &  0.258        \\
				R-VIO  					& 3.502   	& 4.637   & {0.029}   & 0.169          \\
				DSO  					& fail   	& fail   & fail   & fail   \\
				\midrule
			\end{tabular}
			\label{tab:accuracy}
			\vspace{-4mm}
		\end{center}
	\end{table}

	\subsection{Robustness Comparison of Pose Estimation in Smoke Environments }
To evaluate the robustness of our method in challenging environments, we set up smokey environments with a smoke generator.
We compare our method with the famous LOAM \cite{loam} algorithm to evaluate the accuracy and robustness.
The reason for selecting LOAM algorithm is that it can achieve robust pose estimation in light smoke environments while visual odometry will fail in such environments.

    In terms of testing the accuracy of trajectories, we use a pre-built point cloud map as a ground-truth map to present the estimated trajectory. If the estimated trajectory is accurate, it will coincide with the map well. In this way, we can qualitatively evaluate the accuracy of the trajectory.
    In terms of the robustness, we use a smoke machine to adjust the amount of smoke and build light smoke environments (Fig. \ref{fig:first_page}(a)) and heavy smoke environments (Fig. \ref{fig:first_page}(b)) to evaluate the performance of the proposed method. In real scenes, we carried our handheld device and walked from room A to the smoke room by the hallway and returned back to the same position in room A as shown in Fig. \ref{fig:pose_lightsmoke}(a).

	\subsubsection{State estimation comparison under light smoke environments}

	\begin{figure}[htbp]
		\centering
		\includegraphics[width=0.48\textwidth]{experiment_result/odometry_comparison/thermal_vio_hawkins_heavyfog.jpg}
		\caption{The state comparison of our method and LOAM under light smoke environments. (a) is ground-truth map and the red and blue trajectories are generated by LOAM and our algorithm respectively. The green circles show more detailed inspections (a1-a3) on the trajectories and map to illustrate the accuracy. (b) shows the actual point cloud  misalignment from LOAM. (c) shows estimated trajectories in the pre-built map in z direction.}
		\label{fig:pose_lightsmoke}
	\end{figure}

	In this section, we present the trajectory comparison result between our algorithm (blue line) and LOAM \cite{loam} (red line) in light smoke environments.

Fig. \ref{fig:pose_lightsmoke}(a1,a2) shows the trajectories crossing a narrow doorway. Our algorithm (blue trajectory) can pass through the doorway while Lidar trajectories fail and collide with the wall, meaning that Lidar odometry has a relatively large drift. It's worth noting that the narrow doorway has a width of only 0.67 meters.

Fig. \ref{fig:pose_lightsmoke}(a3) shows the result of both trajectories returning to the origin(xyz axes). The red star and yellow star in the figure represent the final position of Lidar odometry and our algorithm, respectively. As we can see, our algorithm's final position is closer to the original start position(xyz axes), indicating that our algorithm is more accurate in this experiment.

Fig. \ref{fig:pose_lightsmoke}(b) shows the actual point cloud misalignment from LOAM under light smoke. This suggests that the Lidar point cloud map has a large drift in such environments. Fig. \ref{fig:pose_lightsmoke}(c) shows the trajectories in the pre-built map in the z-direction, and it can be seen that our algorithm has a small z-directional drift.

	\begin{figure}[htbp]
		\centering
		\includegraphics[width=0.48\textwidth]{experiment_result/odometry_comparison/thermal_vio_heavy_fog.jpg}
		\caption{The state comparison of our method and LOAM under heavy smoke environments. (a) is ground-truth map and the red and blue trajectories are generated by LOAM and our algorithm respectively. The green circles show more detailed inspections (a1-a3) on the trajectories and map to illustrate the accuracy. (b) shows the actual point cloud misalignment from LOAM. (c) shows estimated trajectories in the pre-built map in z-direction.}
		\label{fig:pose_heavysmoke}
	\end{figure}

	\subsubsection{State Estimation Comparison under Heavy Smoke Environments}

	To further test the robustness of our algorithm, we increase the amount of smoke, as shown in Fig. \ref{fig:first_page}(b).

Fig. \ref{fig:pose_heavysmoke}(a1,a2) shows the details of trajectories passing through the narrow doorway, and we can see that our odometry (blue line) can pass through the door while the Lidar odometry (red line) fails and collides with the wall. Fig. \ref{fig:pose_heavysmoke}(a3) shows the result of trajectories returning to the origin (xyz axes), where the red and yellow star represent the final position of Lidar odometry and our algorithm respectively. We can see that the final position of our algorithm is closer to the origin, indicating that our trajectory is more accurate in the heavy smoke environment.
	Fig. \ref{fig:pose_heavysmoke}(b) shows the actual point cloud  misalignment from LOAM \cite{loam} under heavy smoke. This means that Lidar point cloud map has a large drift in such environments. Fig. \ref{fig:pose_heavysmoke}(c) shows the trajectories in the pre-built map in z-direction, and it can be seen that our algorithm has small z-directional drift while Lidar odometry (red trajectory) has large drift. This experiment shows that our algorithm is more robust under heavy smoke.

	\section{Conclusion}
	\label{sec:conclusion}
	In this paper, we propose "ThermalPoint", a CNN-based feature detection network specifically tailored for producing keypoints on thermal images. To achieve real-time performance, we combine ThermalPoint with a novel radiometric feature tracking method and establish reliable correspondences between sequential frames. Finally, a deep feature based thermal-inertial odometry (TP-TIO) estimation framework is proposed and evaluated thoroughly in various visually-degraded environments. Experiments show that our method outperforms state-of-the-art visual and laser odometry methods in smoke-filled environments and achieves competitive accuracy in normal environments.
    In the future, we will add a loop closures module through ThermalPoint and correct accumulated drift for long term navigation.

	\section{ACKNOWLEDGMENT}
This work was partially supported by DARPA agreement \#HR00111820044. Also, this work was supported by National Natural Science Foundation of China (No.61573091), the Fundamental Research Funds for the Central Universities (No.N182608003, No.N172608005), Natural Science Foundation of Liaoning Province (20180520006) and Major Special Science and Technology Project of Liaoning Province (No.2019JH1/10100026).
We also wish to thank Dr. Shehryar Khattak for his constructive advice.

	\balance

	\bibliographystyle{IEEEtran}


\section{Introduction}
The appearance of clusters in light nuclei is well documented~\cite{voe06,Fre18}. The most prominent include the $^{8}$Be ground-state and the $^{12}$C excited Hoyle-state. Both have $\alpha$-particle cluster structures and both play a key role in the synthesis of carbon in stars, through the triple-alpha process. There is good evidence for more complex cluster structures, where, for example, $^{20}$Ne may be described in in terms of $^{16}$O+$\alpha$ substructures and there is a long history of clustering in $^{24}$Mg, where $\alpha$+$^{16}$O+$\alpha$ clustering is one cluster configuration~\cite{voe06}.

The formation of clusters presents an interesting possibility, which is well illustrated by the Beryllium isotopes. The $^8$Be nucleus is unbound to decay to two $\alpha$-particles in its ground-state, whereas $^9$Be is stable. The additional neutron has been found to exist in a delocalised orbit, covalently exchanged between the two $\alpha$-particle clusters. The existence of such structures was discussed by von Oertzen in a series of seminal papers~\cite{voe96,voe97}. These molecular orbitals in the Beryllium isotopes have both $\pi$ and $\sigma$ characteristics, where the notation follows from the atomic analogues. The experimental evidence for such structures comes from the associated rotational bands and decay properties of the states. Perhaps the best example is that of $^{10}$Be, which provides evidence for a molecular structure of two $\alpha$-particles and 2 $\sigma$ covalent neutrons~\cite{Fre06,Suz13}.

The extension of these ideas to three centres and carbon isotopes was made by Milin and von Oertzen~\cite{Mil02}, where the states in $^{13}$C were analysed in terms of their potential molecular characteristics. Convincing evidence for the molecular structure of the carbon isotopes remains to be demonstrated, though there is significant contemporary interest in the molecular structure of $^{14}$C~\cite{Bib19}. Finally, the ideas of nuclear molecules were extended to asymmetric clusters by von Oertzen~\cite{voe01a} and Kimura~\cite{Kim07} through an analysis of the Neon isotopes.

These studies have provided a basis for the understanding of the link between the appearance of molecular orbitals and the underlying cluster structures. This was further cemented through the study of the relationship between the states of the deformed harmonic oscillator and molecular states created through the H\"{u}ckle approach~\cite{McE04}. These molecular structures are found not only in such relatively schematic approaches, but also in more sophisticated methods that treat the nucleons individually, do not impose the clustering a priori and deal with realistic nucleon-nucleon interactions. There are many such approaches, see Ref.~\cite{Fre18} for a review, but the Antisymmeterised Molecular Dynamics, AMD, method has been particularly successful in both reproducing the molecular structures and experimental data, providing a level of confidence that the molecular degrees of freedom are a good representation of the structure of a number of light nuclei.

The present work builds on Ref.~\cite{McE04}, which examined the molecular states of $\alpha$-like nuclei, and extends this understanding to the full spectrum of prolate deformed cluster states across the entire mass range of possible cluster structures.

\section{Theoretical Approach}
The aim of the present paper is to explore the nature of molecular wave-functions across a wide range of nuclear systems. The framework for this is the deformed harmonic oscillator, DHO. The DHO solutions for molecular-like nuclei are  compared with results of calculations where the molecular structures are created explicitly from the orbits of linear combinations of two-centre orbits created using the H\"{u}ckle method.

\begin{figure}
  \begin{center}
  \includegraphics[width=0.38\textwidth]{Fig1a.eps}
    \vspace{-0.0cm}
      \hspace{-0.3cm}
      \caption{Energy levels of the deformed axially symmetric
                 harmonic oscillator as a function
                 of the  quadrupole deformation (oblate and prolate, i.e. negative
                 and positive values of $\epsilon_2$, respectively).
                 Degeneracies (labelled) appear due to crossings of orbits at
                 certain ratios of the length of the long axis (the symmetry axis)
                 to the short perpendicular  axis.
                The regions of high degeneracy define a  shell closure
                 for deformed shapes \cite{Fre95}.}
    \label{fig:HO}
  \end{center}
\end{figure}

\subsection{Deformed Harmonic Oscillator}
The DHO provides a good description of light nuclear systems, where the inclusion of the spin-orbit interaction and an $l^2$ dependence of the nuclear potential results in the Nilsson description of deformed nuclei~\cite{Fre95}. The deformed harmonic oscillator has the well-known shell-like structure where at integer ratios of axial deformations the shell structure is enhanced,~Fig.~\ref{fig:HO}. This shell structure may be associated with stable deformed structures and indeed cluster partitions~\cite{Rae89}.

 The total energy at each level for the DHO is given by:
\begin{equation}
E = \hbar \omega_{x} n_{x} + \hbar \omega_{y} n_{y} + \hbar \omega_{z} n_{z} + \frac{3}{2} \ \hbar \omega_{0},
\end{equation}
where $\omega_{x}$, $\omega_{y}$ and $\omega_{z}$ and $n_{x}$, $n_{y}$ and $n_{z}$ are the oscillation frequencies and number of oscillation quanta in the $x$, $y$, and $z$ directions respectively, and $\omega_{0}$ is the average oscillation frequency. Furthermore, the oscillation frequencies must satisfy:
\begin{equation}
\omega_{0}^3 = \omega_{x}  \omega_{y}  \omega_{z}.
\end{equation}
Therefore, for a nucleus with, for example, a 2:1 deformation along the $z$-direction, the energy levels are given by:
\begin{equation}
E = \frac{6}{5} \ \hbar \omega_{0} n_{x} + \frac{6}{5} \ \hbar \omega_{0} n_{y} + \frac{3}{5} \ \hbar \omega_{0} n_{z} + \frac{3}{2} \  \hbar \omega_{0}.
\end{equation}
The orbitals in the DHO are labelled by [$n_{x}$, $n_{y}$, $n_{z}$], and the number of oscillator quanta in each Cartesian direction determines the orientation of that orbital.

The association of cluster configurations with the shell structure~\cite{Fre95,Rae89} follows from the degeneracies illustrated in Fig.~\ref{fig:HO}. Within the DHO the nucleus $^8$Be, for example, is created through the population of the DHO orbitals,  [$n_{x}$, $n_{y}$, $n_{z}$]$=$[0,0,0]$^4$ and [0,0,1]$^4$.
At a 2:1 deformation this is seen to be associated with the filling of the two orbitals with 2$p$+2$n$, i.e. 2$\alpha$-particles, and hence $^8$Be may be thought of as 2 $\alpha$-clusters. The densities calculated from this DHO configuration also reveal the two-centre cluster structure~\cite{Fre95}. For the next most complex case which is $^{20}$Ne, this also corresponds to a 2:1 shell gap. Here the degeneracies are 2+(2+6) and may be linked to a $^4$He+$^{16}$O cluster structure. The next shell gap at a deformation of 2:1 occurs for degeneracies 2+2+8+8 which would correspond to a 2:1 deformation of $^{32}$S and a $^{16}$O+$^{16}$O cluster structure. A similar situation occurs at a deformation of 3:1, for example the shell gap at cumulative degeneracies of 2+2+2+6+6+6 being linked to the 3:1 deformed state in $^{48}$Cr and a $^{16}$O+$^{16}$O+$^{16}$O cluster configuration.

To form molecular structures in these cluster-systems where the valence particle, typically a neutron, is exchanged between the cluster cores, one needs to examine the properties of the next available orbits beyond those associated with the $\alpha$ and $^{16}$O cluster cores. These are then those which exist immediately above the shell gap. In the case of $^8$Be+$n$ then these are the three levels [1,0,0], [0,1,0] and [0,0,2]. For $^{32}$S it is the 6 orbitals [1,1,0], [2,0,0], [0,2,0], [1,0,2], [0,1,2] and [0,0,4]. In order to ascertain if these have a molecular character, their properties need to be compared with those created explicitly from the available two-centre orbitals. The present work has been framed in terms of neutron covalent orbitals as these are the most prevalent in terms of experimental counterparts. However, equally the ideas may be applied to valence protons or combinations of protons and neutrons.

\subsection{Neutron Wave-functions in Two-Centred Systems}

\paragraph{}The neutron wave-functions may be created from linear combinations of the single-centre obits of the valence neutron around the cluster core. In the case of the $\alpha$-particle the valence orbitals are, [$n_{x}$, $n_{y}$, $n_{z}$]=[1,0,0], [0,1,0] and [0,0,1]. For $^{16}$O, the valence orbitals are the degenerate (in the spherically symmetric HO) [2,0,0], [0,2,0], [0,0,2], [1,1,0], [1,0,1] and [0,1,1] levels. Wave-functions can be linearly combined to create the multiple centre case to compare with the orbital wave-functions found in the DHO. Hence, the nuclear molecular orbit ($\psi$) is then given by a sum of the neutron wave-functions ($\phi$) at each cluster centre, with the appropriate normalisation constants, $c_i$, (determined by the number of cluster centres in the molecule):
\begin{equation}
\psi = c_{1} \phi_{1} + c_{2} \phi_{2} + c_{3} \phi_{3} + \dots
\end{equation}

To find the relative contribution which each single-centre wave-function makes to the overall molecular one,  the H\"{u}ckel method was followed. The appropriate coefficients are found by solving the secular determinant equation which arises from the Schr\"{o}dinger equation for the nuclear molecule (Section~\ref{sec:Huc}).

When making linear combinations of the orbitals in the multiple centre model, the quantum numbers defining the orbital obey certain rules which allowed the combination to be associated with one of the single-centre orbitals~\cite{Fre07, McE04}:
\begin{itemize}
 \item The $n_z$ quantum numbers of the single centre orbits should be summed over the number of centres in the molecule,  $C$, i.e. $n_z=\sum^C_i n_{zi}$, to form the lowest energy molecular configuration.
  \item The $n_x$ and $n_y$ quantum numbers, perpendicular to the direction of deformation, are conserved; the value for $n_{x}$ and $n_{y}$ remains unchanged in the composite system for symmetric clusters ($n_{x1}=n_{x2}=n_{x}$ and $n_{y1}=n_{y2}=n_{y}$).
  \item When forming the orbitals for the  excited molecular states, the quantum number along the direction of deformation can be found using:
  $n_{z} = \sum^{C}_{i} n_{zi} + C - 1$
  where $n_{z}$ is the quantum number for the orbital in the molecule, $n_{zi}$ are the quantum numbers at each centre in the multiple centre molecule.
\end{itemize}

\subsection{The H\"{u}ckel Method}
\label{sec:Huc}
The wave-functions describing all quantum particles can be found by solving the Schr\"{o}dinger equation, for the molecular system
\begin{equation}
\hat{H} \psi = E \psi
\end{equation}
where $\hat{H}$ is the Hamiltonian of the nucleus, $E$ is the energy eigenvalue associated with the particular molecular state, and $\psi$ is the overall wave-function of the nuclear molecule as given above. To demonstrate the method the case of $^{33}$S will be considered. Here there are two $^{16}$O clusters and a single valence neutron. The wave-function describing the valence neutron $^{33}$S is then:
\begin{equation}
\psi = c_{1} \phi_{1} + c_{2} \phi_{2}
\label{eq:wavefunc}
\end{equation}
$\phi_{1}$ and $\phi_{2}$ are the wave-functions of the neutron orbitals at each cluster centre (given by the single centre [2,0,0], [0,2,0], [0,0,2], [1,1,0], [1,0,1] and [0,1,1] orbitals), and the coefficients $c_{1}$ and $c_{2}$ give their relative contributions to the overall wave-function of the nucleus. Equation \ref{eq:wavefunc} is then multiplied by the conjugate wave-function $\psi^{*}$ and integrated over all space, a rearrangement of the equation then gives:
\begin{equation}
E = \frac{\int \psi^{*} \hat{H} \psi \ d\tau}{\int \psi^{*} \psi \ d\tau}
\end{equation}
Substituting the molecular wave-function of $^{33}$S into its eigenvalue equation yields:
\begin{equation}
E = \frac{c^{2}_{1} H_{11} + 2 \ c_{1}c_{2}H_{12} + c^{2}_{2}H_{22}}{c^{2}_1 + 2 \ c_{1}c_{2}S_{12} + c^{2}_{2}}
\end{equation}
Where the equation has been simplified using $H_{ij} = \int \psi_{i}^{*} \hat{H} \psi_{j} \ d\tau$ and $S_{ij} = \int \psi_{i}^{*} \psi_{j} \ d\tau$. The actual wave-function of the nuclear molecule, given by the values of $c_{1}$ and $c_{2}$ will yield the lowest energy state (rather than an excited state); therefore the coefficients can be found by minimising the energy eigenvalue with respect to $c_{1}$ and $c_{2}$. The following secular equations are obtained by finding the partial derivatives and setting them equal to zero:
\begin{equation}
c_{1} \ (H_{11} - E) + c_{2} \ (H_{12} - E S_{12}) = 0,
\end{equation}
\begin{equation}
c_{1} \ (H_{12} - E S_{12}) + c_{2} \ (H_{22} - E) = 0.
\end{equation}
These can be equivalently written in matrix form:
\begin{equation}
\left(
\begin{array}{c c}
H_{11} - E & H_{12} - E S_{12} \\
H_{12} - E S_{12} & H_{22} - E
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\ c_{2}
\end{array}
\right) = \left(\begin{array}{c}
0 \\ 0
\end{array}
\right).
\end{equation}
The energy eigenvalue, $E$, which satisfies this equation is found by solving the secular determinant, because for a non-trivial solution the matrix must be singular:
\begin{equation}
det \left(
\begin{array}{c c}
H_{11} - E & H_{12} - E S_{12} \\
H_{12} - E S_{12} & H_{22} - E
\end{array}
\right) = 0.
\end{equation}

The secular determinant equation may be simplified by making assumptions about the interactions within the nuclear molecule. $H_{11}$ and $H_{22}$ represent the self-interaction energies of the free $^{17}$O nuclei at each centre (i.e. the $n$+$^{16}$O binding energy), which are identical due to the symmetry of the molecule, this value is set equal to $\alpha$. The interaction energy between the two clusters is given by $H_{12} = H_{21} = \beta$ if the $^{16}$O clusters are neighbouring (as is the case for $^{32}$S) or set to zero otherwise. Additionally, due to the neutron wave-functions $\phi_{1}$ and $\phi_{2}$ being located at different cluster centres, their separation enables the simplifying assumption of small overlap to be made, so $S_{12} = S_{21} \simeq 0$. Finally, if the neutron wave-functions at each cluster centre have been normalised then $S_{11} = S_{22} = 1$. Using all of these assumptions gives the secular determinant equation:
\begin{equation}
det \left(
\begin{array}{c c}
\alpha - E & \beta \\
\beta & \alpha - E
\end{array}
\right) = 0.
\end{equation}

The solution to this equation gives $E = \alpha \pm \beta$, and
\begin{equation}
c_{1} \ (\alpha - E) + c_{2} \ \beta = 0 \\
c_{1} \ \beta + c_{2} \ (\alpha - E) = 0,
\end{equation}
the solutions for which are $c_{1}$ = $c_{2}$ and $c_{1}$ = $-\ c_{2}$.

As the molecular wave-function must be normalised, for $^{33}$S the solutions are:
\begin{equation}
\psi = \frac{1}{\sqrt{2}} \ \Big(\phi_{1} \pm \phi_{2}\Big).
\end{equation}
The same method for calculating the molecular wave-function for nuclei comprising more complex arrangements of clusters, i.e. for 3 or more centre systems. The wave-function for each neutron orbital can be constructed in the same way (by solving the secular determinant equation), and linear combinations were made by using the appropriate coefficients. The coefficients for collinear molecules (i.e. when the clusters form a 1d linear structure) of up to six clusters are given in Table \ref{tab1}~\cite{McE04}.

For the two-centre system molecular wave-function is given by:
\begin{equation}
\psi(r) = \frac{1}{\sqrt{2}} \ \Bigg( \frac{1}{\sqrt{\pi^{1/2} 2^{n} n!}}H_{n}(r_{1})exp\left(\frac{- r_{1}^{2}}{2}\right) \pm \frac{1}{\sqrt{\pi^{1/2} 2^{n} n!}}H_{n}(r_{2})exp\left(\frac{- r_{2}^{2}}{2}\right) \Bigg)
\end{equation}
Here, the HO scale factor $\sqrt{m\omega/\hbar}$ is set to unity and $n$ denotes the number of oscillator quanta, $r_{1}$ and $r_{2}$ give the locations of the two clusters with respect to the origin and $H_{n}$ is the appropriate Hermite polynomial. For this investigation, the wave-functions were projected onto the $y-z$ plane and integrated up over the $x$ dimension, and the $z$-direction was chosen to be the deformation axis, or the axis of separation of the clusters.

\begin{table}
\caption{The coefficients of the single particle wave-functions in
the molecular orbits of {\it prolate} nuclei. The magnitude of the
coefficients determine the relative contribution of the
wave-function to the particular molecular orbit, taken from~\cite{McE04}. The increasing orbit number corresponds to higher number of nodes in the composite wave-function and hence higher energy levels.}\label{tab1}
\begin{center}
\begin{tabular} {ccccccc} \hline
 Orbit ($i$) & C$_{i}$ & $^{9}$Be   & $^{13}$C & $^{17}$O & $^{21}$Ne & $^{25}$Mg \\
 \hline
1 & 1 & 0.707           & 0.5       & 0.372   & 0.289 & 0.232    \\
  & 2 & 0.707           & 0.707     & 0.602   & 0.5   & 0.418    \\
  & 3 &                 & 0.5       & 0.602   & 0.577 & 0.521    \\
  & 4 &                 &          & 0.372    & 0.5   & 0.521    \\
  & 5 &                 &          &          & 0.289 & 0.418    \\
  & 6 &                 &          &          &          & 0.232    \\
\hline
2 & 1 & 0.707           & 0.707    & $-$0.602 & 0.5      & $-$0.418 \\
  & 2 & $-$0.707        & 0        & $-$0.372 & 0.5      & $-$0.521 \\
  & 3 &                 & $-$0.707 & 0.372    & 0        & $-$0.232 \\
  & 4 &                 &          & 0.602    & $-$0.5   & 0.232    \\
  & 5 &                 &          &          & $-$0.5   & 0.521    \\
  & 6 &                 &          &          &          & 0.418    \\
\hline
3 & 1 &                 & 0.5   & 0.602    & $-$0.577 & 0.521    \\
  & 2 &                 & $-$0.707    & $-$0.372 & 0        & 0.232    \\
  & 3 &                 & 0.5   & $-$0.372 & 0.577    & $-$0.418 \\
  & 4 &                 &          & 0.602    & 0        & $-$0.418 \\
  & 5 &                 &          &          & $-$0.577 & 0.232    \\
  & 6 &                 &          &          &          & 0.521    \\
\hline
4 & 1 &                 &          & $-$0.372 & 0.5      & $-$0.521 \\
  & 2 &                 &          & 0.602    & $-$0.5   & 0.232    \\
  & 3 &                 &          & $-$0.602 & 0        & 0.418    \\
  & 4 &                 &          & 0.372    & 0.5      & $-$0.418 \\
  & 5 &                 &          &          & $-$0.5   & $-$0.232 \\
  & 6 &                 &          &          &          & 0.521    \\
\hline
5 & 1 &                 &          &          & $-$0.289 & 0.418    \\
  & 2 &                 &          &          & 0.5      & $-$0.521 \\
  & 3 &                 &          &          & $-$0.577 & 0.232    \\
  & 4 &                 &          &          & 0.5      & 0.232    \\
  & 5 &                 &          &          & $-$0.289 & $-$0.521 \\
  & 6 &                 &          &          &          & 0.418    \\
\hline
6 & 1 &                 &          &          &          & $-$0.232 \\
  & 2 &                 &          &          &          & 0.418    \\
  & 3 &                 &          &          &          & $-$0.521 \\
  & 4 &                 &          &          &          & 0.521    \\
  & 5 &                 &          &          &          & $-$0.418 \\
  & 6 &                 &          &          &          &
  0.232 \\
  \hline
\end{tabular}
\end{center}
\end{table}
\section{Results}

\begin{figure}
  \begin{center}
  	\vspace{-4cm}
  \includegraphics[width=0.8\textwidth]{Fig2.pdf}
    \vspace{-4cm}
      \hspace{-0.3cm}
      \caption{\emph{(left)} The deformed harmonic oscillator density calculations for the three nuclei $^{12}$C, $^{16}$O and $^{20}$Ne for the deformations 3:1, 4:1 and 2:1, respectively. \emph{(right)} The calculated densities for the three nuclei $^{12}$C, $^{16}$O and $^{20}$Ne formed from superimposing the densities of the individual clusters, i.e. 3$\alpha$, 4$\alpha$ and $^{16}$O+$\alpha$, respectively.    }
    \label{fig:dho1}
  \end{center}
\end{figure}
\begin{figure}
  \begin{center}
  	\vspace{-4cm}
  \includegraphics[width=0.8\textwidth]{Fig3.pdf}
    \vspace{-4cm}
      \hspace{-0.3cm}
      \caption{\emph{(left)} The deformed harmonic oscillator densities for the nuclei $^{24}$Mg, $^{32}$S and $^{36}$Ar, which reveal the underlying cluster structures $\alpha$+$^{16}$O+$\alpha$, $^{16}$O+$^{16}$O and $^{16}$O+$\alpha$+$^{16}$O, respectively. \emph{(right)} The densities of the three cluster structures $\alpha$+$^{16}$O+$\alpha$, $^{16}$O+$^{16}$O and $^{16}$O+$\alpha$+$^{16}$O created by superimposing the corresponding densities.    }
    \label{fig:dho2}
  \end{center}
\end{figure}
\begin{figure}
  \begin{center}
  	\vspace{-4cm}
  \includegraphics[width=0.8\textwidth]{Fig4.pdf}
   \vspace{-4cm}
      \hspace{-0.3cm}
      \caption{\emph{(left)} The deformed harmonic oscillator densities for the nuclei $^{48}$Cr, $^{56}$Ni and $^{64}$Ge, which reveal the underlying cluster structures $^{16}$O+$^{16}$O+$^{16}$O, $^{40}$Ca+$^{16}$O (asymmetric) and $^{16}$O+$^{16}$O+$^{16}$O+$^{16}$O, respectively. \emph{(right)} The densities of the three cluster structures $^{16}$O+$^{16}$O+$^{16}$O, $^{40}$Ca+$^{16}$O (asymmetric) and $^{16}$O+$^{16}$O+$^{16}$O+$^{16}$O created by superimposing the corresponding densities.}
    \label{fig:dho3}
  \end{center}
\end{figure}
\begin{figure}
  \begin{center}
  	\vspace{-4cm}
  \includegraphics[width=0.8\textwidth]{Fig5.pdf}
   \vspace{-4cm}
      \hspace{-0.3cm}
      \caption{\emph{(left)} The deformed harmonic oscillator densities for the nuclei $^{72}$Kr, $^{80}$Zr and $^{96}$Cd, which reveal the underlying cluster structures $^{16}$O+$^{40}$Ca+$^{16}$O, $^{40}$Ca+$^{40}$Ca and $^{40}$Ca+$^{16}$O+$^{40}$Ca, respectively. \emph{(right)} The densities of the three cluster structures $^{16}$O+$^{40}$Ca+$^{16}$O, $^{40}$Ca+$^{40}$Ca and $^{40}$Ca+$^{16}$O+$^{40}$Ca created by superimposing the corresponding densities.}
    \label{fig:dho4}
  \end{center}
\end{figure}
Figures~\ref{fig:dho1} to ~\ref{fig:dho4} show the densities of the cluster cores built from both the DHO (single-centred nuclear core) and also by explicitly creating a superposition of the densities from the cluster components in a multi-centre approach. These cores correspond to the underlying structures around which the molecular neutrons would orbit. What is clear is that the symmetries that are found in the case of the degeneracies of the DHO (Fig.~\ref{fig:HO}) are broadly reproduced in terms of the densities, revealed in the multi-centre calculations. The multi-centre calculations were performed for separations of the cluster components which are approximately equivalent to the peaks in the density distributions indicated in the DHO densities. The oscillator scale parameter was adjusted from unity for the individual clusters according to the relationship $m\omega/\hbar=0.9A^{1/3}+0.7$ fm$^2$. It can be observed that where the densities of the clusters in the multicentre calculations overlap the densities increase and the details of the multicentre structure are diminished. This can be seen, for example, in the case of the 4$\alpha$ structure of $^{16}$O, in Fig.~\ref{fig:dho1}, where the 4$\alpha$ clusters are most evident in the low density region. However, in the deformed harmonic oscillator the cluster structure is often more evident, demonstrating the importance of the Pauli Exclusion Principle, PEP, in enhancing the cluster structure. Here the overlapping densities from the separated cluster cores does not obey the PEP. This is explored further in the latter part of this paper.

Fig.~\ref{fig:dho1} shows the two sets of calculations for the three isotopes $^{12}$C, $^{16}$O and $^{20}$Ne, which demonstrate the expected cluster structure associated with the 3:1, 4:1 and 2:1 shell closures, respectively. In the case of $^{20}$Ne the multi-centred calculations are for $^{16}$O+$\alpha$  which is reflection antisymmetric.  Similarly Fig.~\ref{fig:dho2} shows the $\alpha$-$^{16}$O-$\alpha$, $^{16}$O+$^{16}$O and $^{16}$O+$\alpha$+$^{16}$O structures in $^{24}$Mg, $^{32}$S and $^{36}$Ar, respectively, associated with the 3:1, 2:1 and 3:1 shell gaps. Evidence for the cluster structure in $^{48}$Cr ($^{16}$O+$^{16}$O+$^{16}$O), $^{56}$Ni ($^{40}$Ca+$^{16}$O), $^{64}$Ge ($^{16}$O+$^{16}$O+$^{16}$O+$^{16}$O), $^{72}$Kr ($^{16}$O+$^{40}$Ca+$^{16}$O), $^{80}$Zr ($^{40}$Ca+$^{40}$Ca) and $^{96}$Cd ($^{40}$Ca+$^{16}$O+$^{40}$Ca) is also found, Fig.~\ref{fig:dho4}, demonstrating that the DHO produces the cluster symmetries. There is no significant reduction of the cluster structure in the deformed harmonic oscillator even in the heaviest systems. As noted, it is these clustered structures which then create the cores for the valence neutrons, covalently exchanged between the clusters.

\subsection{Symmetric Two Centre Molecules}
The $^{9}$Be nucleus has been shown previously to contain two alpha clusters plus a covalent neutron~\cite{McE04}. The molecular characteristics are also found in the AMD calculations of the Beryllium isotopes, described at length in Ref.~\cite{Fre18}. The next two-centre, symmetric, nuclear molecule is $^{33}$S. $^{33}$S can be considered as a 2:1 deformed nucleus with two $^{16}$O clusters with a delocalised neutron. In the DHO model, where the oscillation frequency in the $z$-direction is half that in the $x$ or $y$ directions, the orbitals occupied linked to $^{32}$S are $[n_x,n_y,n_z]=$
[0,0,0], [0,0,1], [0,0,2], [0,1,0], [1,0,0], [0,0,3], [0,1,1] and [1,0,1].
The orbitals available, associated with the next shell closure, for the valence neutron are
[0,0,4], [0,1,2], [1,0,2], [0,2,0], [2,0,0] and [1,1,0], these are all degenerate.
In the two-centre description, each $^{16}$O cluster occupies the orbitals
[0,0,0], [1,0,0], [0,1,0] and [0,0,1],
which leaves the following orbitals available for the valence neutron centred on each cluster $[n_x,n_y,n_z]=$
[0,0,2], [0,2,0], [2,0,0], [0,1,1], [1,0,1] and [1,1,0].

Given that the neutron could exist at either cluster centre with equal probability,  the coefficients associated with the linear combination of the orbitals are $1/\sqrt{2}$ (Eqn. 15). Using the 'rules' described earlier,
the resulting two-centre molecular wave-functions can be matched to the DHO orbitals with a one-to-one correspondence:
\begin{eqnarray}
 \left[0,0,4\right] & \equiv & \frac{1}{\sqrt{2}} \Big[ [0,0,2] + [0,0,2] \Big], \nonumber \\
 \left[0,1,2\right] & \equiv & \frac{1}{\sqrt{2}} \Big[ [0,1,1] + [0,1,1] \Big], \nonumber \\
 \left[1,0,2\right] & \equiv & \frac{1}{\sqrt{2}} \Big[ [1,0,1] + [1,0,1] \Big], \nonumber \\
  \left[0,2,0\right] & \equiv & \frac{1}{\sqrt{2}} \Big[ [0,2,0] + [0,2,0] \Big], \nonumber \\
  \left[2,0,0\right] & \equiv & \frac{1}{\sqrt{2}} \Big[ [2,0,0] + [2,0,0] \Big], \nonumber \\
  \left[1,1,0\right] & \equiv & \frac{1}{\sqrt{2}} \Big[ [1,1,0] + [1,1,0] \Big]. \nonumber \\ \label{eq17}
\end{eqnarray}
Here the $n_z$ quantum number in the DHO is the sum of the two-centre $n_z$ quantum numbers and the $n_x$ and $n_y$ values remain unchanged. All of these linear combinations correspond to the addition (+) of the two-centre orbitals. The negative ($-$) possibility produces the DHO $n_z$ quantum number 2$n_z$+1, i.e. an additional node in the wave-function and hence higher energy orbitals, which can be matched to DHO counterparts in the next oscillator shell at this deformation; but not shown here.
\begin{figure}
  \begin{center}
  	\vspace{-4cm}
  \includegraphics[width=1.\textwidth]{Fig6.pdf}
    \vspace{-6cm}
      \caption{Comparison between 5 of the 6 possible  molecular orbitals for the valence neutron in $^{33}$S calculated in the DHO \emph{(left)} and the multi-centre, H\"{u}ckle, method \emph{(right)}. The [2,0,0] orbital has not been shown. The patterns in the wave-function amplitudes match, demonstrating that the DHO valence orbitals all have a molecular multi-centre nature. The orbitals are labelled by their harmonic oscillator quantum numbers $[n_x,n_y,n_z]$.}
    \label{fig:33S}
  \end{center}
\end{figure}
 Fig.~\ref{fig:33S}  shows the corresponding wave-functions.

 In each case there is a correspondence between the patterns of the wave-function in terms of shape and nodes, though of cause there is not a precise agreement. For example, in the case of the [0,0,4] DHO wave-function, there is good agreement with the molecular wave-function calculated from the linear combination of the two axially displaced [0,0,2] orbitals. Here the final $n_z$ quantum number, 4, is the sum of the two separated $n_z$ quantum numbers, 2+2. This agreement across all molecular wave-functions demonstrates the molecular character of the valence orbitals available for the neutron in $^{32}$S above the shell gap associated with the $^{16}$O+$^{16}$O cluster structure.

The next symmetric two-centred molecule to be examined is $^{81}$Zr, which can be thought of as two $^{40}$Ca clusters with a single delocalised neutron. $^{40}$Ca is the cluster formed in the spherical HO model when orbitals are filled to the next shell closure, after that of the $^{16}$O cluster. In the 2:1 DHO model, the orbitals filled with nucleons from the $^{80}$Zr core are $[n_x,n_y,n_z]=$
[0,0,0], [0,0,1], [0,0,2], [0,1,0], [1,0,0], [0,0,3], [0,1,1], [1,0,1], [0,0,4], [0,1,2], [1,0,2], [0,2,0], [2,0,0], [1,1,0], [0,0,5], [0,1,3], [1,0,3], [0,2,1], [2,0,1] and [1,1,1]. The next highest energy orbitals available for the delocalised neutron are then
[0,0,6], [0,1,4], [1,0,4], [1,1,2], [0,2,2], [2,0,2], [0,3,0], [3,0,0], [1,2,0] and [2,1,0]. These are all degenerate at the 2:1 deformation.
In the two-centre approach, each $^{40}$Ca cluster occupies the orbitals:
[0,0,0], [1,0,0], [0,1,0], [0,0,1], [2,0,0], [0,2,0], [0,0,2], [1,1,0], [1,0,1] and [0,1,1].
The valence neutron then occupies the orbitals
[0,0,3], [0,3,0], [3,0,0], [1,1,2], [1,0,2], [1,1,1], [0,2,1], [2,0,1], [1,2,0] and [2,1,0] around each centre.

As the delocalised neutron in $^{81}$Zr exists at a much higher energy than $^{33}$S, the degeneracy is much larger. This means that there are a greater number of molecular orbitals which must be constructed using linear combinations, but there is still only one possible combination of wave-functions from each cluster centre which give the correct orbital shape, only a subset of these orbitals are shown in Fig.~\ref{fig:81Zr}, but all combinations have been checked. Once again, it can be seen there is a good match between the properties of the 2:1 deformed orbitals and the explicit two-centre molecular counterparts. These two examples demonstrate that at a deformation of 2:1 the deformed harmonic oscillator orbitals above the deformed shell closure all have a molecular structure and thus the deformed harmonic oscillator explicitly includes the 2-centre molecular behaviour.
\begin{figure}
  \begin{center}
  	\vspace{-4cm}
  \includegraphics[width=1\textwidth]{Fig7.pdf}
    \vspace{-6cm}
      \caption{Comparison between selected molecular orbitals (5 of the 10 available possibilities) for the valence neutron in $^{81}$Zr calculated in the DHO \emph{(left)} and the multi-centre, H\"{u}ckle, method \emph{(right)}. The patterns in the wave-function amplitudes match, demonstrating that the DHO valence orbitals have a molecular multi-centre nature. The orbitals are labelled by their harmonic oscillator quantum numbers $[n_x,n_y,n_z]$.}
    \label{fig:81Zr}
  \end{center}
\end{figure}

\section{Three, and More, Centre Nuclear Molecules}
In this section we will use selected examples to show that the behaviour found at 2:1 is reproduced at a deformation of 3:1 and that the inclusion of orbitals with a molecular character is more widely embedded within the deformed harmonic oscillator.

The most simple nuclear molecule with three centres is the collinear arrangement of $^{12}$C, comprising three alpha clusters. The nucleus can also be described with a 3:1 deformation in the DHO model, where the oscillation frequency along the $z$-direction is a third of that along the $x$ or $y$ directions. Following the approach in the previous section, the orbitals filled, by the core nucleons, when the nucleus is considered to exist at a single centre are [0,0,0], [0,0,1] and [0,0,2]. The next available orbitals, to be occupied by the valence neutron in $^{13}$C are
[0,0,3], [0,1,0], [1,0,0].

When the $^{13}$C nucleus is considered as three separate alpha particles, bound by a delocalised neutron, each alpha particle's nucleons occupy orbital [0,0,0]. So, the neutron could be found in a $p$-state above any of the three centres, and correspond to the degenerate
[1,0,0], [0,1,0] and [0,0,1] spherical HO orbits.
\begin{figure}[htb]
  \begin{center}
  	\vspace{-2cm}
  \includegraphics[width=0.5\textwidth]{Fig8.pdf}
\vspace{-3cm}
      \hspace{-0.cm}
      \caption{Amplitudes of two of the three wave-functions for the orbits occupied by valence neutrons in $^{13}$C for \emph{left}) the deformed HO and \emph{right}) the multi-centre orbits.}
    \label{fig:13C}
  \end{center}
\end{figure}
\paragraph{}The coefficients found using the H\"{u}ckel method, when linear combinations of these neutron wave-functions were made, the overall molecular wave-functions possess the form:
\begin{equation}
\psi_1  =  \frac{1}{2}\  \phi_{1} + \frac{1}{\sqrt{2}} \ \phi_{2} + \frac{1}{2} \ \phi_{3} \\
\end{equation}
\begin{equation}
\psi_2  =  \frac{1}{\sqrt{2}} \ \phi_{1} - \frac{1}{\sqrt{2}} \ \phi_{3} \\ \label{eq:mid}
\end{equation}
\begin{equation}
\psi_3  =  \frac{1}{2}\  \phi_{1} - \frac{1}{\sqrt{2}} \ \phi_{2} + \frac{1}{2} \ \phi_{3} \\
\end{equation}

The quantum numbers parallel to the deformation axis, $n_z$, are summed from each centre to give the overall number; the quantum numbers perpendicular to the deformation direction, $n_x$ and $n_y$, were conserved. The results are shown in Fig.~\ref{fig:13C}.

The lowest energy configuration for the [0,0,3] orbital comes from adding the three oscillation quanta situated at each cluster centre in the $z$-direction. Here the phases of the three wave-functions are chosen to match in the overlap region maintaining the number of total nodes in the composite wave-function as 3$\times n_z$.
There is another linear combination in which the central wave-function gives no contribution, Eqn.~\ref{eq:mid}, which produces a wave-function which resembles that of the DHO [0,0,3] orbital. Here an additional node is created in the composite wave-function by the opposite phase of the two wave-functions, generating an additional node and an additional HO oscillator quantum; $2n_z+1$.
For the [0,1,0] orbital the combination of three  oscillator quanta in the $y$-direction yields one overall oscillator quanta perpendicular to the deformation axis. This corresponds to the 3 single-centre wave-functions being added such that the phases combine constructively with no additional nodes in the $z$-direction.

The next three-centre (or 3:1 deformed) nucleus with 3 $^{16}$O cores and a valence neutron is $^{49}$Cr. The possible orbitals in which the valence neutron may reside are
[0,0,6], [0,2,0], [2,0,0], [1,1,0], [0,1,3] and [1,0,3]. Fig.~\ref{fig:Cr49} shows the comparison of 4 of the 6 wave-functions constructed from the DHO and the linear combinations of the neutron wave-function at each of the $^{16}$O clusters. Again, there is a match between the properties of the valence orbitals for the covalent neutron in both the single centred, DHO, and multi-centred, H\"{u}ckle, approaches.

\begin{figure}[htb]
  \begin{center}
  	\vspace{-4cm}
  \includegraphics[width=1\textwidth]{Fig9.pdf}

   \vspace{-6cm}
      \hspace{-0.cm}
      \caption{Amplitudes of the wave-functions for 4 of the 6 possible  orbitals occupied by valence neutrons in $^{49}$Cr for \emph{left}) the deformed HO and \emph{right}) the multi-centre orbitals.}
    \label{fig:Cr49}
  \end{center}
\end{figure}

The structures at a deformation of 4:1 have also been explored in exactly the same way as at 2:1 and 3:1 and the same conclusions emerge; namely there is a 1 to 1 match between the nature of the wave-functions associated with the next available energy levels beyond closing the cluster cores and those predicted in the multi-centred approach. In other words, for \emph{all} structures with identical cluster cores the DHO contains the molecular symmetries for the valence neutrons.

\section{Asymmetric Nuclear Molecules}
Up to this stage, the nuclear molecules considered all have clusters of the same type, e.g. $^{32}$S=$^{16}$O+$^{16}$O and $^{80}$Zr=$^{40}$Ca+$^{40}$Ca clusters. Within the DHO energy level diagram (Fig.~\ref{fig:HO}), this means that the nuclear core will fill orbitals at each centre with the same degeneracies, e.g. 2, 6, 12, ... , and same energies; the valence neutron must then exist in a shell closure which has a new value of degeneracy.
However, nuclear molecules which have a different types of cluster may appear in many systems, e.g. $^{20}$Ne$\equiv^{16}$O+$\alpha$. As pointed out in~\cite{voe06}, it is also possible to construct molecular orbitals from the valence neutron orbital around the two cores, e.g. from $^5$He and $^{17}$O. In this instance, the simplest nuclear molecule of this kind is $^{21}$Ne.

For the nucleus $^{21}$Ne at a deformation of 2:1, the orbitals that are available to the valence neutron are
[0,0,3], [0,1,1] and [1,0,1]. In the multi-centre picture, there are different possible neutron orbitals which could be combined. At one centre the nuclear cluster is an alpha particle, and at the other centre the nuclear cluster is a $^{16}$O nucleus. The different possible orbitals at each centre which should be considered are therfore
[1,0,0], [0,1,0] and [0,0,1], in combination with [2,0,0], [0,2,0], [0,0,2], [1,1,0], [1,0,1] and [0,1,1].

\begin{figure}[htb]
  \begin{center}
  	\vspace{-4cm}
  \includegraphics[width=1\textwidth]{Fig10.pdf}
    \vspace{-6cm}
      \hspace{-0.cm}
      \caption{Amplitudes of  orbitals occupied by valence neutrons in $^{21}$Ne for \emph{left}) the deformed HO and \emph{right}) the multi-centre orbitals, with the $^{16}$O positioned on the left hand side and the $\alpha$-particle on the right.}
    \label{fig:21Ne}
  \end{center}
\end{figure}

\begin{figure}[htb]
  \begin{center}
  	\vspace{-4cm}
  \includegraphics[width=1\textwidth]{fig11.pdf}

   \vspace{-6cm}
      \hspace{-0.cm}
      \caption{Examples of the amplitudes of orbitals occupied by valence neutrons in $^{57}$Ni for \emph{left}) the deformed HO and \emph{right}) the multi-centre orbits. Four cases are shown in the top three panels. }
    \label{fig:57Ni}
  \end{center}
\end{figure}

The linear combinations made to produce the single-centre wave-function follow the same rules in terms of quantum numbers as in the case of identical cluster cores. The quantum number parallel to the direction of deformation ($z$-axis) is found by summing the $n_{z}$ values of the orbitals at each centre, and the quantum numbers perpendicular to the deformation are conserved/unchanged. However, to create wave-functions that have counterparts in the DHO, the transverse quantum numbers, in the $x$ and $y$ directions must match for the multi-centred wave-functions. So for example, as illustrated in Fig.~\ref{fig:21Ne} the linear combination of the [0,1,0] neutron orbit around the $\alpha$-cluster and [0,1,1] neutron orbit about the $^{16}$O-core results in the two centred orbit which has the representation of [0,1,1] within the DHO. If the $\alpha$-particle lies to the right-hand-side of the $^{16}$O then the [0,1,0] wave-function is superimposed on the right-hand component of the clover-like wave-function of the [0,1,1] orbit (see Fig.~\ref{fig:21Ne}).

\begin{table}
	\caption{The relationship between the DHO orbitals for the valence neutron in the $^{40}$Ca+$^{16}$O+n, $^{57}$Ni, system and the orbitals associated with the $^{40}$Ca
		and $^{16}$O clusters.}\label{tab57Ni}
	\begin{center}
		\begin{tabular} {c|cc} \hline
Orbital in $^{57}$Ni & Orbital in $^{40}$Ca & Orbital in $^{16}$O \\
$\left[n_x, n_y, n_z\right]$ & [$n_x$, $n_y$, $n_z$] & [$n_x$, $n_y$, $n_z$] \\
\hline
$\left[0,0,5\right]$ & [0,0,3] & [0,0,2] \\
$\left[1,0,3\right]$ & [1,0,2] & [1,0,1] \\
$\left[0,1,3\right]$ & [0,1,2] & [0,1,1] \\
$\left[2,0,1\right]$ & [2,0,1] & [2,0,0] \\
$\left[0,2,1\right]$ & [0,2,1] & [0,2,0] \\
\hline
		\end{tabular}
	\end{center}
\end{table}

The $^{16}$O+$\alpha$ system is mass-asymmetric and correspondingly the resulting molecular wave-functions do not possess good parity. Parity can be restored by creating linear combinations of the molecular wave-functions formed from having the $\alpha$-particle on either the left or right side of the $^{16}$O core. The resulting wave-functions then more closely resemble the symmetric wave-functions of the deformed harmonic oscillator.

Another asymmetric two-centre nuclear molecule that was investigated was $^{57}$Ni, which could be thought of as  $^{40}$Ca and  $^{16}$O clusters, with a delocalised neutron, at a deformation of 2:1 in the deformed harmonic oscillator. Again, it is possible to construct DHO orbitals from combinations of the single-centre orbitals around each cluster. For example, the valence [0,0,5] DHO orbital in $^{57}$Ni at 2:1 deformation is composed of the [0,0,3] and [0,0,2] orbitals around the $^{40}$Ca and $^{16}$O cores, see Fig.~\ref{fig:57Ni}.  The mapping of the molecular orbitals in $^{57}$Ni to those around the $^{40}$Ca and $^{16}$O clusters in shown in Table \ref{tab57Ni}.

It is also possible to build nuclear molecules from unequal clusters but into symmetric systems. For example, $^{25}$Mg can be constructed from the clusters $\alpha$+$^{16}$O+$\alpha$ with a valence neutron. Such systems were also explored and found that again the wave-functions formed from the linear combination of those at the three centres created the available valence orbital for the neutron in the single centre. Indeed, the ordering of the clusters did not need to be symmetric to achieve this result. Examples explored include  $^{73}$Kr (two $^{16}$O clusters and a $^{40}$Ca cluster) and $^{97}$Cd (two $^{40}$Ca clusters and a $^{16}$O cluster).
\begin{table}
	\caption{The relationship between the DHO orbitals for the valence neutron in the $^{16}$O+$^{40}$Ca+$^{16}$O+n, $^{73}$Kr, system and the orbitals associated with the $^{40}$Ca
		and $^{16}$O clusters.} \label{tab73Kr}
	\begin{center}
		\begin{tabular} {c|ccc} \hline
			Orbital in $^{73}$Kr & Orbital in $^{16}$O & Orbital in $^{40}$Ca & Orbital in $^{16}$O \\
			$\left[n_x, n_y, n_z\right]$ & [$n_x$, $n_y$, $n_z$] & [$n_x$, $n_y$, $n_z$]& [$n_x$, $n_y$, $n_z$] \\
			\hline
			$\left[0,0,7\right]$& [0,0,2] & [0,0,3] & [0,0,2] \\
			$\left[1,0,4\right]$& [1,0,1] & [1,0,2] & [1,0,1] \\
			$\left[0,1,4\right]$& [0,1,1] & [0,1,2] & [0,1,1] \\
			$\left[2,0,1\right]$& [2,0,0] & [2,0,1] & [2,0,0] \\
			$\left[0,2,1\right]$& [0,2,0] & [0,2,1] & [0,2,0] \\
			\hline
		\end{tabular}
	\end{center}
\end{table}
The orbitals for the $^{73}$Kr case are mapped in Table~\ref{tab73Kr}. The result is very similar to that in the case of $^{57}$Ni.

\section{Discussion and Conclusion}

This extensive study of the mapping of molecular orbitals for symmetric and asymmetric cluster systems has found that in every case the molecular orbitals have a one-to-one correspondence with the valence orbitals for neutrons at the deformation given by the underlying cluster structure. In other words, molecular states are not just found for a few select systems but are a ubiquitous feature of nuclear systems where there is an underlying cluster symmetry.

The connection between the deformed harmonic oscillator and the appearance of molecular structures can be traced to the symmetries of the deformed harmonic oscillator~\cite{Naz92} and the ideas of Harvey~\cite{Har75}. The spherical harmonic oscillator possesses a SU(3) symmetry and it has been demonstrated ~\cite{Naz92} that at a deformation of 2:1 and 3:1, etc..., that the deformed harmonic oscillator can be described by 2 or 3 coupled SU(3) groups and hence at a deformation of N:1 the appropriate symmetry is given by N$\otimes$SU(3). In other words, the cluster structures are explicitly encoded in the symmetries. This is seen in the degeneracy patterns in Fig.~\ref{fig:HO} and the densities calculated in Figs.~\ref{fig:dho1} to~\ref{fig:dho4} and Ref.~\cite{Fre95}.
\begin{figure}
	\begin{center}
		\vspace{-0cm}
		\includegraphics[width=0.7\textwidth]{Fig12.pdf}

		\vspace{-2cm}
		\hspace{-0.cm}
		\caption{The Harvey-diagram which demonstrates the connection between the orbitals of the single centre and two-centre system. Infinite separation is when the two clusters are far apart and zero separation is when they are fused.}
		\label{fig:Harv}
	\end{center}
\end{figure}

Harvey~\cite{Har75} noted the rules for combining two oscillator potentials which have been widely used to understand the nature and decay of cluster structures~\cite{Fre07}. These rules relate to the combination of wave-functions from the two centres to the single centre and recognised it is possible to either preserve the total number of nodes in the wave-function or create an additional node at the point where the two wave-functions connect. This means for two wave-functions with harmonic oscillator quantum numbers $n_z$ then two new wave-functions are produced with quantum numbers $2n_z$ and 2$n_z+1$. These rules are depicted in the Harvey Diagram in Fig~\ref{fig:Harv}. For three centres, merging three wave-functions with the same $n_z$ values, produces three possible values of $3n_z$, 3$n_z+1$ and 3$n_z+2$. This can be further extended to N-centres.

For symmetric clusters, the mapping of the clusters observed in the deformed harmonic oscillator to their multi-centre counterparts fully accounts for all of the levels in the separated components and the levels up to the associated shell closure in the deformed harmonic oscillator. The valence orbits above the deformed shell closure map to those outside the closed spherical core.

The merging of the two centres into a single-centre oscillator is a process which respects the Pauli Exclusion Principle, PEP, inasmuch that one of the two merging levels with identical quantum number has an additional node generated. This is also the case for the merging of multi-centres. This effect results in the molecular orbitals not being the low-lying orbitals, which tend to be occupied by the core nucleons, and hence are at higher single-particle energy. The fact that the molecular, and indeed cluster, structures are consistent with the PEP is confirmed through the AMD calculations which explicitly respect the PEP.

For two identical clusters, e.g. two $^{16}$O nuclei, the set of degenerate levels outside the core (e.g. 6 levels in $^{16}$O) map two sets of levels at a deformation of 2:1 split in energy by 3/5$\hbar\omega$ (see Eqn. 3). These set of levels are identical to each other aside from the addition of one unit of $n_z$ (2$n_z$ and 2$n_z$+1). These are the two degeneracies of 12 observed at a deformation of 2:1 in Fig.~\ref{fig:HO}, and are the original two-centre levels which have undergone a 2$n_z$ and 2$n_z$+1 transformation.

Similarly, the set of levels associated with the deformed harmonic cluster structure, up to and including the shell closure, at the deformation of 2:1 are the levels of the separated $^{16}$O clusters that have also undergone the 2$n_z$ and 2$n_z$+1 transformation. So for the original $^{16}$O harmonic orbitals [0,0,0], [0,0,1], [0,1,0] and [1,0,0] these become [0,0,0], [0,0,2], [0,1,0] and [1,0,0] and [0,0,1], [0,0,3], [0,1,1] and [1,0,1]. These are again identical aside from a difference of $n_z$=1. Thus, the 2$n_z$ and 2$n_z$+1 transformation applies equally for the core and valence nucleons. This process of merging from multi-centre to single centre preserves the number of nodes in the final merged wave-function plus creates a new set of wave-functions with an additional node along the separation axis.
These principles apply to more than two centres, e.g. 3 and 4-centre cluster structures, but the transformation will be 3$n_z$, 3$n_z$+1 and 3$n_z$+2.

For asymmetric cores the principle is the same. The merging of the two different valence orbitals around the two different cores preserves the number of total nodes in the final wave-functions and creates a second wave-function with one additional node. These resulting orbitals are the next available orbitals in the composite system, which then have a molecular structure. Linear combinations of orbitals at the two centres should have the same number of transverse nodes to create the molecular wave-functions.

The conclusion from this work is, beyond the observation that all the shell structures in the prolate deformed harmonic oscillator can be associated with cluster structures, that the next available orbitals beyond the shell-gap have molecular structure which is associated with the linear combination of the valence nucleon orbitals around the separated cluster cores. In this way molecular structures should be a wide spread phenomenon and that the molecular behaviour in nuclei is a ubiquitous phenomenon.

The spin-orbit interaction is known to play an important role in the breakdown of clustering, and though it is included in the effective interaction of the AMD approach~\cite{Fre18} where the molecular structures of light nuclei are reproduced, there is no guarantee that the molecular structure persists to much heavier systems described in the present work. The confirmation of the robustness of the states predicted here needs to be tested with more microscopic approaches such as those described in Ref.~\cite{Fre18}.  Moreover, there are experimental challenges for the observation of such states, given that they exist at high excitation energy and will be embedded in a sea of states of other character. This will fragment the molecular states and experimental methods such as that described in Ref.~\cite{Bai19} may be required.

This work is presently being extended to oblate deformations.

\vspace{2cm}