\section{Introduction}
\IEEEPARstart{T}{he} recent success and popularity of Convolutions Neural Networks (CNN) for many computer vision applications have inspired researchers to use those for 3D model classification as well \cite{qi2017pointnet,zhou2018voxelnet,su2015multi}. To exploit the potential of deep networks for this application, different representations of 3D data such as kd-tree \cite{klokov2017escape}, dynamic graphs \cite{wang2019dynamic}, and most recently, spherical harmonics \cite{cohen2018spherical,esteves2018learning,perraudin2019deepsphere,poulenard2019effective} have been proposed. Spherical harmonics is a representation that have attracted significant interest in a wide range of applications including matching and retrieval  \cite{kazhdan2003rotation,wang20163d}, lighting \cite{green2003spherical}, and surface completion \cite{nortje2015spherical}. They attain several favourable characteristics for working with 3D space, such as their basis are defined on the surface of the sphere (volumetric) and are rotation equivariant. In addition, they have shown to provide compact shape descriptors compared to other types of descriptors \cite{kazhdan2003rotation,bulow2001surface}.

The use of CNNs with spherical harmonics has had major success in several recent papers for shape classification \cite{cohen2018spherical,esteves2018learning,poulenard2019effective}, retrieval \cite{esteves2018learning} and alignment \cite{esteves2018learning}. Unlike conventional approaches that use CNNs in regular Euclidean domains, spherical harmonics CNNs (SCNNs) apply convolutions in SO(3) Fourier space, learning features that are SO(3)–equivariant. Interestingly, spherical CNNs have shown to have fewer parameters \cite{cohen2018spherical} and faster training due to the reduction in the dimensionality of the spherical harmonics shape descriptors.

Frameworks that use spherical harmonics CNNs can be divided into two groups: Point-based SCNN that extract features based on point maps or pairwise relations \cite{poulenard2019effective,perraudin2019deepsphere} and the other group that uses spherical harmonics convolution on images casted on the sphere \cite{cohen2018spherical,esteves2018learning}. Approaches that use images casted on the sphere still need to include rotation in the training stage as casted images change with rotation. Nevertheless, it has been shown that spherical convolutional neural networks outperform conventional CNNs, even when both are trained with rotations \cite{cohen2018spherical,esteves2018learning}.

 While spherical CNNs achieved impressive performance in several 3D tasks, their robustness to uncertainty in three-dimensional data has not been explored. Point clouds of 3D models produced by either 3D scanners or multi-view images are often imperfect and contain outliers. As such, development of robust classification frameworks that can deal with such inaccuracies is of interest. Our experiments with using the state-of-art SCNNs with data containing outliers showed that SCNN approaches that are based on point-map features often fail to classify shapes due to outliers distorting the map graph. Also, methods that use image casting typically find the farthest point of intersection between the casted rays and the 3D shape; hence there is a high chance that outlier points are selected as the farthest points. Surprisingly, even when median was used instead of the farthest point, these approaches still often fail to classify shapes.

\begin{figure*}[t]
    \centering
    \includegraphics[scale=.53]{main.pdf}
    \caption{The Proposed Spherical CNN framework. First we sample the shape with $c$ concentric spheres and $n \times n$ grid resolution (seen in orange box). Next we apply Fourier transform (FT) (seen in red box) on the spherical signal and get the basis coefficients up to order \textit{l}. Then we apply the spherical convolution operation where k, is the number of filters. Finally, we feed the spectrum into a fully connected layers (FC) and a classification layer (Cl).}
    \label{fig:fig1}
\end{figure*}

 In this paper, we investigate the robustness of spherical harmonics to data inaccuracies. To start, we explain the spherical harmonics descriptors and the common sampling strategies used in the literature. We then show that using concentric spheres with density occupancy grids provide the highest robustness against data augmentation. The use of concentric spheres generates uniform voxel grids, which have shown to be robust to data augmentation \cite{mukhaimar2019comparative,Riegler_2017_CVPR}. We also propose to use the magnitude of each specific spherical component for shape classification and show that it produces better robustness compared to using the combined magnitudes of different components in each order. In particular, we show that a simple classifier (i.e. fully connected neural network) with the previously mentioned spherical harmonics descriptors and sampling strategy is robust to high levels of data augmentation. Using the above knowledge and the inspiration from the recent success of spherical CNNs approaches ~\cite{cohen2018spherical,esteves2018learning}, we propose a robust spherical convolutional neural network framework (called RSCNN) that is able to deal with different types of uncertainty inherent in three-dimensional data measurement. The framework is outlined in Figure 1. In our approach, the entire spherical convolution operation is performed in the Fourier domain (we show that using the inverse transform reduces robustness to data augmentation), and the sampling is performed over concentric spheres.

Our proposed approach uses voxel grids, similar to 3D CNNs~\cite{Riegler_2017_CVPR}, that have shown to be robust to the influence of outliers. However, we are using spherical grid that has been shown to exhibit better rotation equivarience compared to euclidean grids \cite{boomsma2017spherical} and using spherical CNNs that have less trainable parameters \cite{cohen2018spherical}. Unlike previous approaches \cite{qi2017pointnet,qi2017pointnet++,su2015multi}, the applied spherical convolution operation is simply multiplying the filter kernels by the spherical harmonic coefficients, hence, the convolution operation does not disrupt the input signal through the use of a pooling operation or grid altering. Also, we will show the output features produced by the convolution operation are highly robust.

Our key contributions in this paper are as follows:
\begin{itemize}
\item We investigate the robustness of spherical CNNs and spherical harmonics as it has not discussed before.
\item We propose a new spherical CNN framework that is significantly more robust than previous spherical CNNs approaches.
\item We demonstrate the efficiency and accuracy of our method on shape classification with the presence of several types of data inaccuracies, and we show that our framework significantly outperform all previous approaches.
\end{itemize}

The rest of this paper is structured as follows. An in-depth literature review is provided in section \ref{sec:Related work}, while preliminaries required to explain spherical harmonics are provided in section \ref{sec:Preliminaries}. A detailed discussion of the proposed model is presented in section \ref{sec:Method}. Description of an extensive set of experiments is included in section \ref{sec:Experiments} as well as a comparison with other the state of the art 3D shape classification (including other spherical harmonics) methods. We also examine the robustness of different spherical harmonics structures and descriptors in section \ref{sec:ABLATION STUDY}. Section \ref{sec:Conclusion} concludes the paper.

\section{Related work}
\label{sec:Related work}
Spherical harmonics have been used for 3D shape classification for many years. Early classification frameworks used spherical harmonic coefficients as shape descriptors \cite{kazhdan2003rotation}. Later, with the use of spherical convolutions, classifier networks were empowered to learn descriptive features of objects. We study both approaches in terms of their performance under data augmentation.

The spherical CNNs proposed in \cite{cohen2018spherical,esteves2018learning}, for spherical signals defined on the surface of a sphere, addresses the rotation equavarience using convolutions on the SO3 rotation group. In \cite{cohen2018spherical}, the input spherical signal is convolved with S2 convolution to produce feature maps on SO3, followed by an SO3 convolution,
while in \cite{esteves2018learning}, spherical convolutions were used. Steerable filters  \cite{worrall2017harmonic,weiler2018learning,weiler20183d} were used to achieve rotation equivariance, where filters were restricted to the from of complex circular harmonics \cite{worrall2017harmonic}, or complex valued steerable kernels were used \cite{weiler2018learning}. In \cite{poulenard2019effective}, a network called Sphnet is designed to apply spherical convolution on volumetric functions \cite{atzmon2018point} generated using extension operators applied on point cloud data. Unlike previous approaches, spherical convolution is applied on point clouds instead of spherical voxel grid, resulting in better rotation equivariant. In another work, a network called Deepsphere \cite{perraudin2019deepsphere},  spherical CNNs are used on graph represented shapes, where shapes are projected onto the sphere using HEALPix sampling, in which the relations between the pixels of the sphere build the graph. The graph is then represented by the Laplacian equation, which is solved using spherical CNNs. Ramasinghe et al. \cite{Ramasinghe} investigated the use of radial component in spherical convolutions instead of using spherical convolutions on the sphere surface. They proposed a volumetric convolution operation that was derived from Zernike polynomials. Their results show that the use of volumetric convolution provides better performance by capturing depth features inside the unit ball. Spherical signals have also been used in conjunction with conventional CNNs by \cite{boomsma2017spherical,su2017learning,coors2018spherenet} to achieve better rotational equivariance than signals on euclidean space, but unlike spherical CNNs, their CNN operations are not equivariant. You et al. \cite{you2018prin} used concentricity to sample 3D models with a sampling strategy that has better robustness to rotation.
While previous approaches used spherical harmonics to build rotation equavarient neural networks, our focus is on building a robust spherical harmonics structure. As such, our choice of representation is not restricted to spherical harmonics descriptors that are rotation equavarient.

In terms of recent deep learning approaches for 3D shape classification, 3D CNNs have been used with voxel based 3D models \cite{maturana2015voxnet,wu20153d,zhou2018voxelnet,Riegler_2017_CVPR} using several occupancy grids \cite{maturana2015voxnet}. Such a representation has shown to be robust to data augmentation \cite{Riegler_2017_CVPR} while some implementations (e.g.~\cite{brock2016generative}) have achieved very high classification accuracy for clean objects (using ModelNet40).
Another approach is to use 2D CNN on images of the 3D mesh/CAD objects rendered from different orientations \cite{su2015multi,johns2016pairwise,wang2019dominant}. The rendered images are usually fed into separate 2D CNN layers, a pooling layer follow these layers to aggregate their information. These methods take advantage of existing pertained models to achieve high classification accuracy. When testing MVCNN \cite{su2015multi}, the classification accuracy was heavily affected by data augmentation, especially outliers. Another approach uses unsorted and unprocessed point clouds directly as an input to the network layers \cite{qi2017pointnet,qi2017pointnet++}. These approaches use a max-pooling layer that was tested to be robust to point dropout and noise \cite{mukhaimar2019comparative}. However, when tested with outliers, their performance was significantly affected. Another approach is to build upon relations between points \cite{klokov2017escape}; for such methods, the existence of outliers completely changes the distance graph and causes such an approach to fail.

In terms of robust classification frameworks that exist in the literature, Pl-net3D \cite{mukhaimar2019pl} decompose shapes into planar segments and classify objects based on the segments information. DDN~\cite{gould2019deep} proposes an end-to-end learnable layer that enables optimization techniques to be implemented in conventional deep learning frameworks. An m-estimators based robust pooling was proposed instead of max pooling used in conventional CNNs. Our approach shows better robustness to data argumentation while involves less computation.

\section{Preliminaries}
\label{sec:Preliminaries}

In this section, we review the theory of spherical harmonics along with their associated descriptors that are used for classification tasks. In addition, we review the theory of convolution operations applied to spherical harmonics.

\subsection{Spherical harmonics}

Spherical harmonics are a complete set of orthonormal basis functions that are defined on the surface of unit sphere $S2$ as:

\begin{equation}
Y_{l}^{m}(\theta,\varphi)=\sqrt{{(2l+1)(l-m)!\over 4\pi(l+m)!}}P_{l}^{m}(\cos\theta)e^{im\varphi}
\label{equ:sh}
\end{equation}

where $P_{lm}(x)$  is the associated Legendre polynomial, $l$ is the degree of the frequency and $m$ is the order of every frequency $(l\geq 0, m\leq\vert l\vert). ~\theta\in[0,\pi],   ~\varphi\in[0,2\pi]$ denote the latitude and longitude, respectively. Any spherical function $f(\theta,\varphi)$ defined on unit sphere $S2$ can be estimated by the linear combination of these basis functions:
\begin{equation}
  f(\theta,\varphi) = \sum\limits_{l=0}^{\infty}\sum\limits_{m=-l}^{\l}\hat{f}_{lm}Y_{l}^{m}(\theta,\varphi)
  \label{equ:fun}
\end{equation}

where $\hat{f}_{lm}$ denotes the Fourier coefficient found from:
 \begin{equation}
  \hat{f}(l,m) = \int_{S^{2}}f(\theta,\varphi)\bar{Y_{l}}^{m}(\theta, \varphi)d\varpi
  \label{equ:coef}
\end{equation}

The following descriptors implies that any spherical function can be described in terms of the amount of energy $|\hat{f}|$ it contains at every frequency:
 \begin{equation}
  F1 =(|\hat{f}_{0,0}|,|\hat{f}_{1,0}|,|\hat{f}_{1,1}|,...|\hat{f}_{lm}|)
  \label{equ:desc1}
\end{equation}

 \begin{equation}
  F2 =(|\hat{f}_{0}|,|\hat{f}_{1}|,|\hat{f}_{2}|,...|\hat{f}_{l}|),~where~  \hat{f}_{l} =\sqrt{\sum\limits_{m=-l}^{\l}\hat{f}_{lm}^2}
  \label{equ:desc2}
\end{equation}

Both of these descriptor vectors have been used for classification of shapes \cite{kazhdan2003rotation,wang20163d}. However, only the second descriptor is rotation equavareint while the first one carries more shape information. We have investigated the use of both descriptors for shape classification and the result is provided in section~\ref{sec:Method}.

\subsection{Spherical convolution}

If we have a function $f$ with its Fourier coefficients $\hat{f}$ found from equation~(\ref{equ:coef}), and another function or kernel h with its Fourier coefficients $\hat{h}$, then the convolution operation in spherical harmonic domain is equal to the multiplication of both functions Fourier coefficients as shown below \cite{driscoll1994computing}:

 \begin{equation}
 (f * h)_{l}^{m} =  \sqrt {{4\pi} \over {2\ell +1}} \hat{f}_{\ell}^{m} \hat{h} _{\ell}^{0}.
   \label{equ:fh}
\end{equation}

Here, the convolution at degree $l$ and order $m$ is obtained by multiplying of the coefficient $\hat{f}_{\ell}^{m}$ with the zonal filter kernel $\hat{h}_{\ell}^{0}$. The inverse transform is also achieved by summing over all l values:

\begin{equation}
  (f * h)(\theta,\varphi) = \sum\limits_{l=0}^{\infty}\sum\limits_{\vert m\vert \leq l}(f * h)_{l}^{m}Y_{l}^{m}(\theta,\varphi).
  \label{equ:fh1}
\end{equation}

\section{Method}
\label{sec:Method}
To facilitate the description of our solution space, we first introduce the shape classification problem and the effect of outliers, followed by depiction of our proposed solution architecture (as shown in Figure~\ref{fig:fig1}). The proposed solution starts from sampling the input spherical signal, then moves to convolution operations, and finally performs classification.

\subsection{Problem statement}

The point clouds of 3D models produced by either 3D scanners or multi view images are often imperfect and contain outliers. To be able to mitigate the effect of outliers, we need to develop ways to model the data. Our primary objective for using spherical harmonics is to achieve robust classification rather than rotation equavarience. As such, we constrain the object rotations to a rotation around z-axis similar to most contemporary approaches including \cite{qi2017pointnet,zhou2018voxelnet,wang2019dynamic}.

\subsection{Sampling on concentric spheres}

The first step in using the spherical harmonics for modelling an object is to sample the input signal, which is referred to as $f(\phi,\theta)$ in equation~(\ref{equ:coef}). Two types of sampling are used in literature \cite{kazhdan2003rotation}: Sampling over the concentric spheres, and sampling over the sphere surface. Our results show that the method \cite{esteves2018learning} that use image casting fails to mitigate the effect of outliers. Therefore, the second sampling strategy is adopted in this paper. The generated grids for both types of sampling change with rotation and therefore, rotation needs to be included during training.

We generate a spherical voxel grid that consists of $c$ concentric spheres with $n \times n$ grid resolution for each concentric sphere. The generated spherical voxel grid allows sampling over the unit ball $(S3)$, with each voxel being represented by $(r,\theta,\phi)$ where ($r \,\to\, 0:c.$ $ \theta,\phi \,\to\, 0:n)$. We distribute the given 3D shape over the grid and we keep record of number of points inside each voxel and produce a density occupancy grid. The use of such occupancy grid is expected to provide reliable estimates in the presence of outliers. We compare occupancy grids in the next section.

\subsection{Classification with Fourier coefficients}

The robustness of the two descriptors represented by equation~(\ref{equ:desc1}) and equation~(\ref{equ:desc2}) for object classification are compared in this section by feeding each of them to a fully connected neural network. We tested their robustness against: Gaussian noise with 0.10 standard deviations, uniformly scattered outliers with 50\% percentage, and 80\% Random point dropout. More about those tests are detailed in the experiment section. The results are shown in table~\ref{Tab:t1}. The results in table~\ref{Tab:t1} also show the robustness of the sampling strategies and the occupancy grids mentioned in the previous section. The results show that the used sampling and the density occupancy grid provides a high degree of robustness to outliers. In addition, these results also show that the descriptor in equation~(\ref{equ:desc1}) F1 provides higher classification accuracy than using the descriptor in equation~(\ref{equ:desc2}) F2 as the first one carries more shape information.

\begin{table}[]
 \centering
\begin{threeparttable}
\caption{Classification accuracy results for objects augmented with noise, missing points, and outliers.}
\label{Tab:t1}
 \fontsize{10}{12}\selectfont
\begin{tabular}{|l|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\begin{tabular}[c]{@{}c@{}}occupancy   \\ grid\end{tabular}} & clean  & \begin{tabular}[c]{@{}c@{}} random \\ dropout\end{tabular} & \begin{tabular}[c]{@{}c@{}} noise\end{tabular} & \begin{tabular}[c]{@{}c@{}} outliers\end{tabular} \\ \hline \hline
binary                                                                            & \textbf{0.79} & 0.34                                                           & 0.24                                              & 0.14                                                 \\ \hline
density+F1                                                                        & 0.78   & \textbf{0.75}                                                           & \textbf{0.37}                                            & \textbf{0.50}                                                     \\ \hline
density+F2                                                                        & 0.68  & 0.68                                                         & 0.3                                                 & 0.24                                                 \\ \hline
\end{tabular}
    \begin{tablenotes}
      \small
      \item The dataset details are shown in section~\ref{sec:Experiments}.
    \end{tablenotes}
  \end{threeparttable}

\end{table}

\subsection{Implementation of Spherical convolution}

We propose to apply the spherical convolution on 3D models that are decomposed into concentric spheres. The use of concentric spheres generates a uniform spherical voxel grid that enables the spherical convolution neural network to learn features over the unit ball (as appose to only learning over the unit sphere). We use separate convolution operation at each concentric sphere to allow our network to learn features relevant to that sphere. To achieve the spherical convolution, we use equation~(\ref{equ:fh}) in which the learned kernel is a zonal ($m=0$) filter h with dimension of $l*c$ , where $l$ is the frequency degree, and $c$ is the number of concentric spheres. Similar to \cite{esteves2018learning}, we parameterize the kernel filters in the spectral domain. No inverse Fourier transform is applied after the spherical convolution. Therefore, our convolution operation is entirely in the spectral domain, which reduces the convolution computation time. Our results show that applying inverse Fourier transform (IFT) diminishes the robustness to outliers as the overall accuracy reduces by almost 15 percent. This can be related to equation~(\ref{equ:fh1}), where for each $\theta$ and $\phi$, the output signal is calculated by summing the entire coefficients. Thus, if the coefficients are already altered by outliers, the output signal error will be magnified/accumulated due to this summation. A detail discussion on this topic is provided in the Appendix.

Our experiments show that applying the convolution operation works well with augmented data and the network has been able to learn better features and be more discerning in terms of object classification. This is shown by applying t-sne~\cite{maaten2008visualizing} to clean and augmented data and the results are provided in the Appendix. As the application of convolutions on 3D voxel grids has shown to be robust to the influence of outliers \cite{mukhaimar2019comparative}, we would expect our method to exhibit a high degree of robustness to outliers as well.

Compared to previous approaches, unlike other networks such as pointnet~\cite{qi2017pointnet} where their max pooling chooses outliers as max, our proposed method does not use pooling or grid altering operations. Our proposed convolution operation can be described as follows: Let $x \in X$ be our input spherical coefficients at a given degree, the spherical convolution operation in equation~(\ref{equ:fh}) is simply $f(x_{i})=k \times x_{i}$, where $k$ is a constant calculated from the square root term of the same equation and the kernel value at that degree followed by the non-linearity operation. This mathematical operation does not alter the input signal and only assist with extracting better features in both clean and augmented data as seen from the t-sne results (provided in Appendix B).

Compared to 3D CNNs such as the octnet~\cite{Riegler_2017_CVPR}, spherical CNNs are rotation equavarient, which could help in increasing our performance. In addition, the use of spherical convolution have shown to have less trainable parameters, where one layer is enough to achieve good performance \cite{Ramasinghe}.

\subsection{classification layer}

The returned feature map by the convolution operations represents the feature vector defined in equation 4. The map is then fed to a fully connected and a classification layers.  The feature vector in equation 5 could be used as well, however, it is less robust to data augmentation as shown in table 1. Although the feature vector defined in equation 4 is not rotation invariant, given that we are training with rotations, we would expect our network to learn rotations.

\section{Experiments} \label{sec:Experiments}

The following experiments are carried out with $z$ rotation, only. We compare our framework with the state of the art published spherical convolution architectures, point cloud classification methods, and robust methods. We considered outliers, noise, and missing points as our types of data augmentation in this paper since corruption of point clouds with such inaccuracies are commonplace.

\subsection{Dataset}

To test different methods, we use the benchmark dataset ModelNet40 \cite{qi2017pointnet++}. ModelNet40 dataset contains point clouds of 40 different shapes, with 9,843 point clouds for training, and 2468 for testing. Each object consists of 2048 points, which are normalized within the unit cube. We generated three artificial datasets from the ModelNet40 dataset for each of the three types of data augmentation. These datasets are only used during testing, while we only add noise and random rotations as our data augmentation during training.

\subsection{Architecture}

For the following experiments, we used one convolution layer having 16 output channels. We use prelu after the convolution operation as our non-linearity. The number of concentric spheres was chosen to be 7 with grid resolution of 64 by 64 of each sphere. The basis degree is selected to be 9. We compare different architectures and study the effect of changing the number of spheres in our ablation study.

\subsection{Training}

Similar to \cite{cohen2018spherical,esteves2018learning}, we also include rotation in the training data. The patch size is set to 16, learning rate varies from 0.001 to 0.00004, and number of epochs is set to 48. We used a TITAN Xp GPU, where only 350Mb of memory was used during training.\\

\subsection{Robustness to outliers}

Similar to \cite{mukhaimar2019pl}, we present two outlier scenarios generated with different mechanisms. In the first scenario, we test our model in the presence of scattered outliers: points uniformly distributed in the unit cube. In the second scenario, added outliers are grouped into clusters of ten or twenty points, which are uniformly distributed in the unit cube (similar to \cite{ning2018efficient}). The overall number of scattered points for this scenario are fixed to ten or twenty percent as shown in table~\ref{Tab:t4}. Points in each cluster are normally distributed with zero mean and standard deviations of 4\% and 6\%.

\begin{figure}[]
\centering
\includegraphics[scale=.33]{outl.pdf}

\caption{Classification accuracy versus outliers.}
\label{fig_outl}
\end{figure}

Figure~\ref{fig_outl} and table~\ref{Tab:t4} show that our model is highly robust to the influences of outlier in both scenarios. The classification accuracy only drops by 6\% percent when half the data are outliers, while PL-net3D and other models drop by significantly higher margins. For Spherical-cnn~\cite{esteves2018learning}, even when we used the median aggregation for generating the unit sphere grid (instead of max aggregation), the network remains sensitive to the influences of the outliers. Similarly, SpH-net~\cite{poulenard2019effective} performs poorly when there were outliers as these outliers distort the distance graph. Results for using Spherical-cnn with our occupancy grid are examined in our ablation study.

\begin{table}[h]
\caption{Object classification results on clustered outliers}
\label{Tab:t4}
\centering
 \fontsize{10}{12}\selectfont
\begin{tabular}{lllllllll}
\cline{1-4}
\multicolumn{1}{|c|}{\multirow{3}{*}{Method}} & \multicolumn{1}{c|}{10\%}    & \multicolumn{1}{c|}{10\%}    & \multicolumn{1}{c|}{20\%}    &  &  &  &  &  \\
\multicolumn{1}{|c|}{}                        & \multicolumn{1}{c|}{10p}     & \multicolumn{1}{c|}{10p}     & \multicolumn{1}{c|}{20p}     &  &  &  &  &  \\
\multicolumn{1}{|c|}{}                        & \multicolumn{1}{c|}{$\mathcal{N}$(0.04)} & \multicolumn{1}{c|}{$\mathcal{N}$(0.06)} & \multicolumn{1}{c|}{$\mathcal{N}$(0.04)} &  &  &  &  &  \\ \cline{1-4}
\multicolumn{1}{|c|}{Octnet}                  & \multicolumn{1}{c|}{0.47}    & \multicolumn{1}{c|}{0.48}    & \multicolumn{1}{c|}{0.37}    &  &  &  &  &  \\ \cline{1-4}
\multicolumn{1}{|c|}{PL-Net3D}                    & \multicolumn{1}{c|}{0.79}    & \multicolumn{1}{c|}{0.8}    & \multicolumn{1}{c|}{0.67}    &  &  &  &  &  \\

\cline{1-4}
\multicolumn{1}{|c|}{Ours}                    & \multicolumn{1}{c|}{\textbf{0.81}}    & \multicolumn{1}{c|}{\textbf{0.81}}    & \multicolumn{1}{c|}{\textbf{0.75}}    &  &  &  &  &  \\

\cline{1-4}
\multicolumn{4}{p{170pt}}{ 10\%: outliers percentage, 10p: 10 Points in each cluster.}
\end{tabular}

\end{table}

\subsection{Robustness to noise}

We simulated the effect of noise in point cloud data by perturbing points with zero mean, normally distributed values with standard deviations ranging from .02 to .10 as shown in Figure~\ref{fig_noise}. Compared to other models, our proposed model performance deteriorated the least (by around 18\%) for relatively large amount of noise (at 0.10 noise level). SPHnet was significantly affected by noise, while spherical-CNN performance was relatively much better than SPHnet.

\begin{figure}[]
\centering
\includegraphics[scale=.33]{noise.pdf}

\caption{Classification accuracy versus noise.}
\label{fig_noise}
\end{figure}

\subsection{Robustness to missing points}

We also evaluated the effect of missing points on the classification performance of different methods. Figure~\ref{fig_missp} shows that our model classification accuracy drops by only 2\% when half the points are eliminated, and by 20\% when 90\% of points are removed. spherical-cnn classification accuracy drops by 10\% when half the points are eliminated and it degrades after that. SPHnet classification accuracy drops by 8\% when 60\% of points are eliminated and it degrades after that. Ocntnet classification accuracy degrades after 50\%.

\begin{figure}[]
\centering
\includegraphics[scale=.33]{miss.pdf}

\caption{Classification accuracy versus missing points.}
\label{fig_missp}
\end{figure}

\section{Ablation study} \label{sec:ABLATION STUDY}

We conducted a number of experiments to explore all possible solutions of our method and their performance under different data augmentation and the results are shown in table~\ref{Tab:comp}. The first row shows the results of our proposed model (RSCNN) with 3D Models having 2000 points (same results shown in the previous section), while the second and the rest of the raw's show the results when training without noise addition. We implemented the inverse transform operation after applying the convolution in our model. As a result, our model performed worse and became less robust to data augmentation.  The outputs are presented in the third row of table~\ref{Tab:comp}. The effect of performing Inverse Fourier Transform on the classification accuracy is discussed in the Appendix.
In the next step, we evaluated our model with no fully connected layer to reduce the number of trainable parameters. However, the results shown in the third row, suggest that such an action is detrimental for the overall performance.

As most existing techniques used only 2000 points in their experiments, we also used the same number of points in our comparative experiments. However, for objects with 10,000 points, as shown in the forth row, our network achieves higher performance and robustness. The classification accuracy with 50\% outliers increased to 79\%, while we noticed 1\% enhancement for both clean and noisy objects.
When testing Spherical-cnn~\cite{esteves2018learning} utilizing the same grid we used, the performance of their method enhanced tremendously. Although spherical-cnn was able to achieve 70\% classification accuracy with 50\% outliers, it is still around 10 percent lower than the performance of the proposed method.

\begin{table}[]
\caption{Classification accuracy versus different network architectures and different data augmentations}
\label{Tab:comp}
\fontsize{10}{12}\selectfont

\begin{tabular}{|l|c|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{method}} & \multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}\\  sampled\\  points\end{tabular}} & \multicolumn{4}{c|}{classification accuracy}                                                                                \\ \cline{3-6}
\multicolumn{1}{|c|}{}                        &                                                                                & clean & \begin{tabular}[c]{@{}c@{}}80\% \\  dropout\end{tabular} & \begin{tabular}[c]{@{}c@{}}0.1\\ noise\end{tabular} & \begin{tabular}[c]{@{}c@{}}50\%\\ outliers\end{tabular} \\ \hline
ours                                          & 2000                                                                                & \textbf{0.82}  & 0.72                                                           & 0.62                                                & 0.74                                                    \\ \hline

ours*                                        & 2000                                                                                & 0.81  & 0.72                                                           & 0.63                                                & 0.75                                                    \\ \hline

ours + IFT                                    & 2000                                                                                & 0.78  & 0.71                                                           & 0.45                                                & 0.58                                                    \\ \hline
ours no FC                                    & 2000                                                                                & 0.75  & 0.65                                                           & 0.54                                                & 0.46                                                    \\ \hline
\multicolumn{6}{|l|}{}                                                                                \\[-2ex] \hline
ours                                          & 10000                                                                                & \textbf{0.82}  & \textbf{0.80}                                                            & \textbf{0.64}                                                & \textbf{0.79}                                                    \\ \hline
sph-cnn \cite{esteves2018learning}*                                & 10000                                                                                & 0.80   & 0.53                                                           & 0.43                                                & 0.69                                                    \\ \hline
\end{tabular}
\end{table}

\subsection{Grid resolution}

We tested our method with 4, 5, 7, and 10 concentric spheres where each sphere has a 64 by 64 grids. The results are shown in Figure~\ref{fig_c}. For both clean and outlier-augmented objects, our network performance gradually increases up to using 7 concentric spheres and plateaus afterwards. As such, we used only 7 concentric spheres in our experiments.

\begin{figure}
\centering
\includegraphics[scale=.3]{grid1.pdf}

\caption{Classification accuracy versus number of concentric spheres.}
\label{fig_c}
\end{figure}

\subsection{Number of layers}

We tested with different architectures having one to four layers as shown in Figure~\ref{fig_l}. We find that increasing the number of layers from 1 to 4 does not improves the accuracy on clean and augmented objects, thus we restricted our model to a single convolutional layer.

\begin{figure}[]
\centering
\includegraphics[scale=.33]{num_layers.pdf}

\caption{Classification accuracy versus number of layers for our network.}
\label{fig_l}
\end{figure}

\section{Limitation and future work}

 In theory, we expect our model to achieve better performance compared to Spherical-cnn~\cite{esteves2018learning} as we learn features from the unit ball. However, Spherical-cnn~\cite{esteves2018learning} was able to achieve higher classification accuracy for clean objects. Nevertheless, The proposed framework performance is competitive to other robust methods that exist in the literature \cite{mukhaimar2019pl,gould2019deep}.

\section{Conclusion}  \label{sec:Conclusion}

Spherical harmonics have been used for object classification over the years with several frameworks existing in the literature. Their rotation equavarience, and the reduced number of trainable parameters are key properties they exhibit. However, their robustness to data augmentation have not been studied before.
Point clouds of 3D models are often imperfect and contain outliers. Therefore, it is important to use a robust classification framework that can overcome such inaccuracies.
In this paper we present the robustness of spherical harmonics approaches available in the literature for classification of point cloud objects. We also propose a robust model that have shown to compete with the state-of-art models. Our model uses the voxel grid of concentric spheres to learn features over the unit ball. In addition, we keep the convolution operations in the Fourier domain without applying the inverse transform used in previous approaches. As a result, our model is able to learned features that are less sensitive to data augmentation. We tested our proposed model against several types of data augmentation such as noise and outliers. Our results show that the proposed model outperforms the state of art networks in terms of robustness to effects of data augmentation.

\appendices
\section{Effect of inverse Fourier transform }

Our experiments show that applying inverse Fourier transform (IFT) reduces robustness to outliers by almost 15 percent. To demonstrate this, we compared the convolution operation of both cases in Figure~\ref{fig_case}, where~(a) represents the convolution operation that uses IFT, and~(b) represents our convolution operation. In Figure~\ref{fig_case}-(a), we plot the input voxel grid as a 2D image having a dimension of 64 by 64*7, where 7 is the number of concentric spheres. The pixel values represent the values of the voxel grids. The color bar shows how many points inside each voxel.  While in Figure~\ref{fig_case}-(b), we plot the Fourier coefficients instead of the voxel grid. For each case, a signal of a clean shape and a signal of the same shape having outliers are used as inputs. Next we apply the convolution and plot their outputs in the second column. Since in (a) the convolution uses IFT, the output signal is again returned to the Euclidian domain. However for our convolution in (b), our output signal is kept in the Fourier domain.
Finally, we plot the difference between both outputs in the third column with its associated histogram in the last column. Comparing the histograms of (a) and (b) show that our convolution operation produces less error.

\begin{figure*}[]
\centering
\includegraphics[scale=.7]{case_study.pdf}

\caption{A comparison between the spherical convolution operation with (a) IFT applied after the convolution, and (b) no IFT applied. In both scenarios, we have a clean shape and the same shape augmented with outliers.  The first column is the input signal, where in (a) we have a voxel grid signal that consists of 64*64*7 voxels (spherical grid). The color bar represents the number of points inside each voxel. while the first column in (b) represents the spherical harmonics coefficients of the same input signal shown in (a). The second column shows the output of the spherical convolution operation, where in (a) we apply IFT, hence, the signal is in the spherical grid domain, while in (b) the signal is kept in the Fourier domain. The third column shows the difference in signal between the clean and outlier case, and the last column shows the histograms of the images in the third column. Histograms have the same bin width.}
\label{fig_case}
\end{figure*}

\section{t-SNE}
We used the t-Distributed Stochastic Neighbor Embedding (t-
SNE) \cite{maaten2008visualizing} to visualize 2D embedding of learnt shape features in figure~\ref{fig_tsne}. We only show a few representative labels results. The t-sne technique helps us to visualizes high-dimensional feature vectors (1024 dimensions) aggregated from the convolution operation into the 2D space shown in the figure. Figure~\ref{fig_tsne}-(a,b) show the t-sne embedding for a simple classifier that uses equation~(\ref{equ:desc1}) for classification.  We show scenarios of clean data of ModelNet40 test split, and the 50\% outlier augmented data respectively. Similarly in figure~\ref{fig_tsne}-(c,d) for our model.
The results show that using our model helped in better grouping similar labels for clean data and augmented data as well. We still can see the labels clustered clearly in the presence of outliers, which means that our learned feature is almost not affected by the presence of outliers.

\begin{figure*}[]
\centering
\includegraphics[scale=.7]{tsne.pdf}

\caption{Visualization of 2D embedding of shape features using t-SNE technique for: (a) a simple classifier using clean data, (b) same classifier  using 50\% outliers augmented data, (c) our model on clean data, and (d) our model on 50\% outliers augmented data. }
\label{fig_tsne}
\end{figure*}

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}


\section{\label{sec:Introduction}Introduction}

In-plane structural anisotropy in individual layers of van der Waals materials produces electronic bandstructures that are unique in contrast to the family of Transition Metal Dichalcogenides (TMDCs) with in-plane rotational symmetry. Rhenium (Re) based  Group VII TMDCs like ReS\textsubscript{2} have garnered considerable attention because of their intriguing anisotropic optical, vibrational and electronic properties arising from reduced crystal symmetry \cite{aslanLinearlyPolarizedExcitons2016a,tongayMonolayerBehaviourBulk2014, aroraHighlyAnisotropicInPlane2017a,echeverryTheoreticalInvestigationsAnisotropic2018,jadczakExcitonBindingEnergy2019b,urbanNonEquilibriumAnisotropic2018,simUltrafastQuantumBeats2018a,olivaPressureDependenceDirect2019,liuIntegratedDigitalInverters2015a,cuiTransientAbsorptionMeasurements2015,zhouStackingOrderDriven2020,wangDirectIndirectExciton2019,simSelectivelyTunableOptical2016,hafeezRheniumDichalcogenidesReX2017}. This originates from ReS\textsubscript{2}’s distorted 1T structure \cite{wilsonTransitionMetalDichalcogenides1969}, where the extra electron from the Re atom  contributes to the strong Re-Re metal bond, forming  a zigzag chain along the b-axis. This results in higher electron mobility along the b-axis \cite{liuIntegratedDigitalInverters2015a}, and optical anisotropy manifesting itself as highly polarized photoluminescence (PL) and absorption due to two strongly bound exciton species X\textsubscript{1} and X\textsubscript{2}  with dipole moments along different in-plane directions \cite{aslanLinearlyPolarizedExcitons2016a,jadczakExcitonBindingEnergy2019b,urbanNonEquilibriumAnisotropic2018,simUltrafastQuantumBeats2018a}. These properties open the door to a class of polarization-sensitive, on-chip devices like  polarization controlled all-optical switches \cite{simLightPolarizationControlledConversion2019}, polarized LEDs \cite{wangPolarizedLightEmitting2020}, photodetectors \cite{zhangReS2BasedFieldEffectTransistors2015} and polarization-based quantum logic gates \cite{kwonAll2DReS2Transistors2019}. Additionally, for Group VI TMDCs the electronic and optical properties are strongly dependent on the number of layers of the crystal, and show a drastic change in the monolayer limit \cite{splendianiEmergingPhotoluminescenceMonolayer2010}. Conversely, ReS\textsubscript{2} with its weak interlayer coupling shows no such drastic change when going from bulk to monolayer \cite{tongayMonolayerBehaviourBulk2014}, making it ideal for multilayer photonic device applications. In order to create highly sensitive devices, however, a meticulous study of its anisotropic dielectric properties is required. Furthermore, there has been considerable debate about the nature of ReS\textsubscript{2}’s bandgap, with the recent consensus being that it is marginally indirect except in bilayer form \cite{urbanNonEquilibriumAnisotropic2018,olivaPressureDependenceDirect2019,wangDirectIndirectExciton2019,biswasNarrowbandAnisotropicElectronic2017,webbElectronicBandStructure2017}. Experimental studies that can shed light into this matter are essential.

\begin{figure*}[t]
	\includegraphics{fig1}
	\caption{\label{fig:1}(a) Top and side view of the distorted 1T phase of monolayer  ReS\textsubscript{2}. The arrow indicates the b-axis of crystal orientation. b) Optical microscope image of exfoliated  ReS\textsubscript{2} transferred on Si/SiO\textsubscript{2} substrate. The white scale bar represents 10 $\mu$m. (left). Optical contrast along the dotted yellow line showing layer thickness (right). (c) Differential reflectance from unpolarized reflectance study at 3 different points in the sample marked as 1, 2 and 3 (2, 4 and 11 nm) in (b). Further studies are conducted on point 3, i.e., 11 nm ReS\textsubscript{2}. (d) Schematic showing the exciton complexes at Z and K1 points of the Brillouin Zone along with the low-lying momentum-dark states (e) Differential reflectance for two different incident polarizations w.r.t the b-axis as indicated in the inset of the top and bottom panel, showing four exciton peaks for 11 nm ReS\textsubscript{2}. (f) Polarized PL intensity measured at two different orientations of the analyzer w.r.t the b-axis as indicated in the inset of the top and bottom panel.  (e), (f) corresponds to the polarization state for which the contribution is minimized from X\textsubscript{2} and X\textsubscript{4} (top panel) and X\textsubscript{1} and X\textsubscript{3} (bottom panel) respectively. X\textsubscript{R} denotes the Rydberg excitations of X\textsubscript{1} and X\textsubscript{2}.}
\end{figure*}

In this work, we focus on high-resolution, polarization-resolved reflectance and PL at low temperature, discovering two additional exciton peaks that were predicted for ReS\textsubscript{2} \cite{echeverryTheoreticalInvestigationsAnisotropic2018} but not observed in previous studies. It is envisaged that these two shoulder peaks denoted by X\textsubscript{3} and X\textsubscript{4} are observed due to the splitting of spin degenerate excitonic states by combined effect of electron-hole exchange interaction, structural anisotropy, and spin-orbit coupling.  X\textsubscript{3} and X\textsubscript{4} appear on the higher energy sides in both the absorption and PL measurements with similar polarization orientation as X\textsubscript{1} and X\textsubscript{2} respectively. The anisotropic nature of ReS\textsubscript{2}’s dielectric properties is further demonstrated via polarization dependent reflectance. Transfer matrix method was utilized to understand the asymmetric excitonic lineshape in our reflectance from the ReS\textsubscript{2}/SiO\textsubscript{2}/Si dielectric stack as a function of layer thickness and thus extract the real and imaginary part of anisotropic refractive index of the material. In addition, we report excitation polarization dependence of the PL, which provides insight into the anisotropic absorption of this material. For excitation energy 1.88 eV, the integrated PL intensity from all four exciton species is maximized when excitation polarization is oriented along X\textsubscript{2}. Furthermore, our temperature-dependent study finds the four exciton peak positions are well resolved up to a temperature of 150 K, and their dipole orientations are preserved. Most intriguingly, we observe an anomalous PL intensity variation around 30 K, which necessitated the development of a model using rate equations, considering phonon scattering from low-lying dark exciton states. Our model provides evidence of the existence of indirect exciton states which are closely-lying below the bright excitons, confirming the quasi-indirect nature of ReS\textsubscript{2}’s bandgap. Our analysis can shed light on the understanding of the excitonic properties in this anisotropic material and provide useful information for photonics and optoelectronic device engineering with ReS\textsubscript{2}.

\section{\label{sec:Results}Results and Discussion}

\begin{figure*}[t]
	\includegraphics{fig2}
	\caption{\label{fig:2}(a) Colour plot representing the differential reflectance spectrum as a function of incident polarization angle $\theta_{P}$ measured w.r.t. the b-axis for 11 nm ReS\textsubscript{2}. Line profile in the lower panel shows reflectance at a particular polarization angle where all four excitons can be observed. (b) Polar plots of differential reflectance for four different excitons X\textsubscript{1}, X\textsubscript{2}, X\textsubscript{3} and X\textsubscript{4} proportional to the oscillator strength (values for X\textsubscript{3} and X\textsubscript{4} are multiplied by 3 for clarity). (c) Polarization resolved absolute reflectance data (black) along with the simulation result (red) obtained from the transfer matrix method, shown for three different polarisation angles. (d), (e) Colour plots demonstrating the in-plane anisotropy in the real and imaginary part of the refractive index, obtained from the fitting of absolute reflectance at variable polarization. The lower panel shows the corresponding line profile obtained at particular angle of polarization.}
\end{figure*}

A schematic top view is shown in Fig. \ref{fig:1}(a) of a single layer of ReS\textsubscript{2}, where each molecular layer is a sandwich of a Re layer between two S layers. The direction of the b-axis is marked by a solid arrow which can be identified from the optical image during measurement. Sample was prepared by mechanical exfoliation technique, and dry transferred on to 340 nm SiO\textsubscript{2}/Si substrate. The optical image of the sample is shown in Fig. 1(b) where we can approximately identify the number of layers from the optical contrast along a line profile, which has been further verified using AFM (see Fig. S1 in the Supplemental Material \cite{SeeSupplementalMaterial}). On-chip crystallographic orientation can be estimated by observing the sample edges, utilizing the fact that ReS\textsubscript{2} is likely to cleave along the axis containing the covalent Re-S bonds, which is parallel to the b axis \cite{aslanLinearlyPolarizedExcitons2016a,liangOpticalAnisotropyAudoped2009a}.

Reflectance measurement was performed using a broadband halogen source with a spot size of $\sim$3 \textmu m at three different points of the sample as shown in Fig. 1(b): Points 1, 2 and 3, being 2 nm, 4 nm and 11 nm thick respectively. At 11 nm thickness (point 3), ReS\textsubscript{2} approaches its bulk character, wherein PL intensity becomes independent of the number of layers \cite{aslanLinearlyPolarizedExcitons2016a,tongayMonolayerBehaviourBulk2014}. However, its 2D character is preserved since excitons are confined mostly in single layers \cite{aroraHighlyAnisotropicInPlane2017a}.

\begin{figure*}[t]
	\includegraphics{fig3}
	\caption{\label{fig:3}(a) Colour plot showing the PL spectrum in the log scale as a function of excitation laser polarization angle $\theta_{P}$ w.r.t. the b-axis for 11 nm ReS\textsubscript{2} when no analyzer is placed at the output path. The dotted line indicates the polarization angle at which the overall PL intensity is maximum. The linear scale PL spectrum when the excitation laser aligned with the b-axis is shown in the lower panel.  (b) Colour plot showing the PL spectrum in the log scale as a function of angle of the analyzer (w.r.t. b-axis) obtained at a fixed polarization of the excitation laser which maximize the PL intensity. For the same excitation polarization, a linear scale PL spectrum at 30$^\circ$ angle of the analyser is shown in the lower panel. (c) Polar plot representing the PL intensity variation for two most prominent peaks, X\textsubscript{1} and X\textsubscript{2} as a function of the excitation polarization angle. (d) Polar plot representing the PL intensity variation of X\textsubscript{1}, X\textsubscript{2}, X\textsubscript{3}, and X\textsubscript{4} as a function of the angle of the analyzer $\theta_{P}$ at fixed excitation polarization (Intensity of X\textsubscript{1}, X\textsubscript{3} and X\textsubscript{4} are multiplied by factors of 2, 3 and 3 respectively for clarity). The dotted lines indicate corresponding angle of maximum PL for X\textsubscript{1} and X\textsubscript{2} excitons.}
\end{figure*}

Excitonic features in PL and reflectivity are too broad at room temperature to resolve all the peaks, therefore, the sample was cooled down to 3.2 K in a closed cycle cryostat. The reflectivity measured at point 1 and 2 show an asymmetric lineshape about the excitonic resonance, as shown in Fig. 1(c). We see a pronounced red shift of exciton peak positions from few layers to bulk limit in agreement with previous reports \cite{aslanLinearlyPolarizedExcitons2016a,jadczakExcitonBindingEnergy2019b,urbanNonEquilibriumAnisotropic2018}. In reflectivity data, we plot the differential reflectance, that is, $\frac{\Delta R}{R}=(R_{ReS_2+SiO_2/Si}-R_{SiO2/Si})/R_{SiO2/Si}$, where $R_{ReS_2+SiO_2/Si}$ and $R_{SiO_2/Si}$ are the reflectance spectra from the sample and the SiO\textsubscript{2}/Si substrate respectively. Due to interference in the multilayer film system consisting of Si, SiO\textsubscript{2} and ReS\textsubscript{2} as shown in the schematic in Fig. 1(a), a broad anti-reflection dip coincides with the excitonic resonances. This makes the 11 nm ReS\textsubscript{2} (at point 3) ideal to probe well-resolved exciton peaks in reflectance, which almost mimics the PL spectrum. Chosen thickness of the SiO\textsubscript{2} layer and the dielectric properties of 11 nm ReS\textsubscript{2} make these peaks show a background free, almost symmetric Fano lineshape. In contrast, the reflectivity measured at point 1 and 2 show asymmetric Fano lineshape, which is a result of interference with the background reflectivity with a relative phase change of $\pi$ across the exciton resonance \cite{limonovFanoResonancesPhotonics2017}. The Fano asymmetry parameter varies across the three points. The peaks at higher energy range are the higher-order Rydberg series of the excitons \cite{aslanLinearlyPolarizedExcitons2016a,jadczakExcitonBindingEnergy2019b}. What follows below are based on experimental results conducted on the 11 nm ReS\textsubscript{2}.

We observe four peaks at 1.545$\pm$0.001 (X\textsubscript{1}), 1.558$\pm$0.001 (X\textsubscript{3}), 1.576$\pm$0.001 (X\textsubscript{2}) and 1.593$\pm$0.001 (X\textsubscript{4}) eV. The peaks at higher energy range are the higher-order Rydberg series of the excitons, denoted by X\textsubscript{2} \cite{aslanLinearlyPolarizedExcitons2016a,jadczakExcitonBindingEnergy2019b}. The newly observed peaks X\textsubscript{3} and X\textsubscript{4} are in contrast to what has been observed earlier in ReSe2 \cite{aroraHighlyAnisotropicInPlane2017a}, since they appear at the higher energy side of X\textsubscript{1} and X\textsubscript{2} respectively.  However, from the resemblance between these additional peaks and the similar peaks for ReSe2, we speculate their origin is the splitting of singlet and triplet states of excitons due to electron-hole exchange interaction. Other plausible contributions for such splitting are the broken rotational symmetry due to structural anisotropy and spin-orbit coupling in ReS\textsubscript{2}.  As shown in Fig. 1(d), we attribute the higher energy prominent exciton transition X\textsubscript{2} to the K1 point, and X\textsubscript{1} to the Z point of the Brillouin Zone, using the results of a recent, comprehensive ab-initio calculation \cite{olivaPressureDependenceDirect2019}. X\textsubscript{D1} and X\textsubscript{D2} are low-lying momentum-dark states nearly degenerate with X\textsubscript{1} and X\textsubscript{2} respectively. They will be of interest when considering temperature dependence of PL intensity later in this paper. Figures 1(e) and 1(f) show the polarization resolved differential reflectance and PL at a polarization angle of -10$^{\circ}$  (top panels) and 90$^{\circ}$  (bottom panels) respectively, where all exciton peaks as seen in the reflectivity are also observed in PL with same energy positions. The 2 meV Stokes shift between PL emission and absorption peak (Fig. S2 in the Supplemental Material \cite{SeeSupplementalMaterial}) and the strong PL intensity indicates the pristine quality of our sample. We see no change in its optical properties over time and over multiple cooling cycles (Fig. S3 \cite{SeeSupplementalMaterial}).

In the polarization resolved reflectance measurement, linearly polarized white light was used as the source and the reflected light was collected via 0.7 NA objective lens, while the excitation polarization angle $\theta_{P}$ is varied with respect to (w.r.t.) the b-axis using a half-wave plate. The experimental result for 11 nm ReS\textsubscript{2} is shown in the colour plot in Fig. \ref{fig:2}(a). A line plot in the bottom panel is chosen at an angle where we see all the four exciton resonances. The polar plot Fig. 2(b) shows that the two exciton species X\textsubscript{1} and X\textsubscript{2} are polarized at angles of 6$^{\circ}$ and 79$^{\circ}$  respectively w.r.t the b-axis (For fitting function see Supplemental Material \cite{SeeSupplementalMaterial}), which agrees with earlier reported values within experimental error \cite{aslanLinearlyPolarizedExcitons2016a,olivaPressureDependenceDirect2019,hoCompleteseriesExcitonicDipole2019}. X\textsubscript{3} and X\textsubscript{4} follow the polarization variation of X\textsubscript{1} and X\textsubscript{2} respectively.

As shown in Fig. 1(a), we model the polarization resolved reflectance via transfer matrix technique \cite{blakeMakingGrapheneVisible2007,liMeasurementOpticalDielectric2014a,zhangMeasuringRefractiveIndex2015,hsuThicknessDependentRefractive2019,kuzmenkoKramersKronigConstrained2005}  (see Supplemental Material \cite{SeeSupplementalMaterial}) to obtain the reflectance spectrum of the ReS\textsubscript{2}/SiO\textsubscript{2}/Si stack and hence extract the frequency dependent dielectric function of ReS\textsubscript{2}. We fit the polarization-resolved reflectance to obtain the dielectric function $\epsilon(\omega)$ for every polarization. The dielectric function is given by $\epsilon(\omega)=\epsilon_b+\sum\frac{f_i}{\omega_{0i}^2-\omega^2-i\omega\gamma_i}$, where $\epsilon_b$, $f_i$, $\omega_{0i}$ and $\gamma_i$ are the background dielectric constant, oscillator strength, resonance frequency and the linewidth of the i\textsuperscript{th} oscillator, and the summation is over all exciton resonances.  Fig. 2(c) shows absolute reflectance spectra obtained at three different polarization angles of the linearly polarized incident beam, along with their theoretical fits. All the fitted parameters are given in Tables S1-S4 in Supplemental Material  \cite{SeeSupplementalMaterial}.

The anisotropic refractive index $\tilde{n}\left(\hbar\omega\right)=n\left(\hbar\omega\right)+ik\left(\hbar\omega\right)$, where, $n$ and $k$ are the real and imaginary parts of refractive index, are shown in the colour plots Figs. 2(d) and 2(e) respectively. Interestingly, as shown in Fig. 2(c), we observed that the transfer matrix model best fit our experimental data with the polarization resolved reflectance only. However, in the Supplemental Material (Fig. S4 \cite{SeeSupplementalMaterial}) we show that it is not possible to obtain a good fit for unpolarized reflectance data with this model. As the unpolarized reflectance is a result of averaging over all polarization, it cannot be fitted by an effective value of the dielectric constant, since it is in truth highly anisotropic.

A 660 nm laser with a spot size of $\sim$1 \textmu m is used to excite the sample with variable polarization w.r.t the b-axis. At first, the integrated PL was collected directly at the spectrometer slit without any analyzer at the output port. We discover that the overall PL intensity varies with the excitation polarization direction as shown in the colour plot in Fig. \ref{fig:3}(a). The intensity of the overall spectrum is modulated; however, the intensity ratio of X\textsubscript{1} and X\textsubscript{2} remains unchanged as the incident polarization is varied. Intensity from all exciton species is maximum at a particular polarization of excitation laser, corresponding to an in-plane direction which turns out to be the polarization direction of the X\textsubscript{2} exciton, within experimental error. This is a consequence of absorption at excitation energy 1.88 eV also being anisotropic, peaking at the angle along X\textsubscript{2}. This can be predicted from the angle dependent absorption obtained from our fitting (Fig. S5 \cite{SeeSupplementalMaterial}), and agrees with theoretical calculations \cite{echeverryTheoreticalInvestigationsAnisotropic2018}. The incident polarization dependence for two other lower excitation energies (1.76 and 1.70 eV) closer to X\textsubscript{2} resonance was also tested, which showed slightly different behavior (Fig. S6 \cite{SeeSupplementalMaterial}). The anisotropy becomes more pronounced as the excitation energy comes closer to the exciton resonances.

\begin{figure*}[t]
	\includegraphics{fig4}
	\caption{\label{fig:4}(a) PL spectra measured at different temperatures for a fixed angle of the analyser and excitation polarization. The dotted lines are a guide to eye. (b) Schematic diagram for the rate equation model. Curved arrows indicate all the phonon scattering processes considered in the model. Black dotted lines indicate non-radiative decay channel. Yellow region indicates the light cone for momentum allowed optical transitions. The solid straight arrows represent radiative decay channels. (c) Variation of PL intensity with temperature showing a local maxima around 30 K. Blue and red solid lines are fits obtained from the rate equation model. (d) Peak positions for the four excitons as a function of temperature. Solid lines are the fitted curves from the theoretical model.}
\end{figure*}

Next, we keep the excitation polarization fixed at an angle for which the overall PL is maximum and measure the polarization state of the PL via an analyzer placed before the spectrometer slit. Data shown in Fig. 3(b) was recorded while rotating the analyzer angle $\theta_{A}$ from 0$^{\circ}$  to {360}$^{\circ}$, where zero is along the b-axis. From the polar plot in Fig. 3(d), we find the two dominant emission peaks at X\textsubscript{1} and X\textsubscript{2} are polarized w.r.t the b-axis at angles of 6$^{\circ}$ and {82}$^{\circ}$ respectively, which is same as the angles obtained from reflectance considering experimental error. It is observed as shown in Figs. 2(b) and 3(d) that the two shoulder peaks at X\textsubscript{3} and X\textsubscript{4} follow the same polarization directions as X\textsubscript{1} and X\textsubscript{2} respectively.

We perform a temperature-dependent PL measurement, the results of which are plotted in Fig. \ref{fig:4}(a). The analyzer is fixed throughout at an angle for which all four exciton peaks are visible. The two shoulder peaks X\textsubscript{3} and X\textsubscript{4} are not resolvable above 150 K since they are dominated by X\textsubscript{1} and X\textsubscript{2} whose linewidths increase with temperature. We observe an anomalous temperature variation of PL intensity as shown in Fig. 4(c), with a maxima in intensity near 30K, which can be attributed to the presence of low-lying dark states (X\textsubscript{D1}, X\textsubscript{D2}). Spin-forbidden dark excitons  \cite{molasBrighteningDarkExcitons2017} have recently been reported to be the cause of anomalous  temperature dependent PL in Group VI TMDCs  \cite{aroraDarkTrionsGovern2020}. On the other hand, momentum-forbidden low-lying states have been hypothesized before to explain the ratio of X\textsubscript{1} and X\textsubscript{2} populations not following the expected Boltzmann distribution \cite{aslanLinearlyPolarizedExcitons2016a,urbanNonEquilibriumAnisotropic2018}. This is also evident in Fig. \ref{fig:3} where it is observed that the X\textsubscript{1} PL intensity is lower than X\textsubscript{2} which is persistent even at higher temperatures as shown in Fig. 4(c).  Time-resolved measurements have also revealed that the radiative lifetimes of excitons to be less than 10 ps \cite{wangDirectIndirectExciton2019}, indicating the excitons are not thermalized.

To gain a comprehensive understanding of the temperature variation of this hot photoluminescence, we propose a rate equation model involving the bright (X\textsubscript{1}, X\textsubscript{2}) and dark states (X\textsubscript{D1}, X\textsubscript{D2}) as shown in the schematic diagram Fig. 4(b), where excitons in these states are generated via continuous wave (CW) pumping. Two kinds of phonon scattering are important in this model which are discussed below. First, scattering of excitons from X\textsubscript{D1} to X\textsubscript{1} and X\textsubscript{D2} to X\textsubscript{2} with rates $\gamma_1$ and $\gamma_2$ respectively, where $E_{01}$ and $E_{02}$ are the energy of the phonons involved. Second, $\alpha_1$, $\alpha_2$ and $\beta_1$, $\beta_2$ are the scattering rates from X\textsubscript{1} and X\textsubscript{2} to momentum-forbidden dark states outside the light cone. Excitons scatter to a continuum of states outside the light cone via phonons of all available energies. To simplify the model, we have considered only two different phonon energies per exciton state - $E_{1n}$, $E_{1m}$ for X\textsubscript{1} and $E_{2n}$, $E_{2m}$ for X\textsubscript{2}. This effectively accounts for the scattering at both low temperature and high temperature regimes (see Fig. S10 in the Supplemental Material \cite{SeeSupplementalMaterial}).

We assume $\tau_1$ and $\tau_2$ are the radiative recombination times from X\textsubscript{1} and X\textsubscript{2}, and $\tau_{D1}$ and $\tau_{D2}$  are the non-radiative recombination times from X\textsubscript{D1} and X\textsubscript{D2}. By solving the rate equation under CW excitation (see Supplemental Material \cite{SeeSupplementalMaterial}), the PL intensity for X\textsubscript{1} and X\textsubscript{2} is obtained as a function of temperature, and fitted with the experimental data as shown in Fig. 4(c). From this model we find the two momentum-dark states X\textsubscript{D1} and X\textsubscript{D2} are present 16 meV and 5 meV below X\textsubscript{1} and X\textsubscript{2} states respectively. These dark states strongly indicate a quasi-indirect band gap at Z and K1 point in the Brillouin zone. This is in agreement with ab-initio calculations which indicate ReS\textsubscript{2} is marginally indirect \cite{urbanNonEquilibriumAnisotropic2018,olivaPressureDependenceDirect2019}.

When temperature is initially increased (3.2 K$<$T$<$30 K) X\textsubscript{1} and X\textsubscript{2} are thermally populated from X\textsubscript{D1} and X\textsubscript{D2} states by absorbing phonons. At the same time, X\textsubscript{1} and X\textsubscript{2} excitons are scattered from the radiative window to outside the light cone, where they recombine non-radiatively, by absorbing phonons of energy 17 meV and 3.2 meV respectively. In this temperature regime, we observe dark states are appreciably more populated than bright states on laser excitation, which results in a net gain in excitons scattered to the bright states. This causes PL intensity to increase to a maximum at 30K. At higher temperatures, the scattering process is mainly dominated by phonons of energies 52 meV and 41 meV from X\textsubscript{1} and X\textsubscript{2} states respectively. Scattering to dark states increases more rapidly than that to bright states after 30 K, causing PL intensity to decrease. All phonon energies obtained, except for 5 and 3.2 meV, have recently been reported for ReS\textsubscript{2} \cite{mccrearyIntricateResonantRaman2017a}. From our fitting we also notice that the effective scattering time from X\textsubscript{1} and X\textsubscript{2} is much more than their recombination times, even up to high temperatures, which corroborates the fact that X\textsubscript{1} and X\textsubscript{2} emit hot photoluminescence.

We plot the observed redshift of exciton peak positions with temperature, and fit with the model describing the temperature dependence of a semiconductor bandgap \cite{odonnellTemperatureDependenceSemiconductor1991}: $E_g\left(T\right)=E_g\left(0\right)-S\hbar\omega\left[\coth{\left(\frac{\hbar\omega}{2kT}\right)}-1\right]$. Here $E_g\left(0\right)$ is the exciton resonance energy at T = 0 K, S is a dimensionless electron-phonon coupling constant, and $\hbar\omega$ is the average phonon energy. The fitted parameters are $E_g\left(0\right)$, $S$ and $\hbar\omega$. The values of $S$ are 1.74, 1.78, 2.42 and 1.83 for X\textsubscript{1}, X\textsubscript{2}, X\textsubscript{3} and X\textsubscript{4} respectively, and $\hbar\omega$ is around 20$\pm$3 meV, concurring with earlier reports \cite{hoAbsorptionedgeAnisotropyReS1998,hoOpticalStudyStructural2005}. Energy separation between the four exciton peaks remains constant with temperature as can be seen in Fig. 4(d). The angle between the excitons’ dipole moment orientation and b-axis remains unchanged with temperature indicates the temperature independent anisotropy (Fig. S7 \cite{SeeSupplementalMaterial}). Finally, temperature dependent reflectance measurements were used to acquire the corresponding dielectric function using the transfer matrix method (Fig. S8 \cite{SeeSupplementalMaterial}). $\epsilon_1$ shows significant change only near the exciton resonances.

\section{\label{sec:Conclusion}Conclusion}

In conclusion, we observed two additional exciton shoulder peaks at higher energy sides of X\textsubscript{1} and X\textsubscript{2} which are attributed to splitting of spin degenerate exciton states. The ab initio calculations required to systematically investigate the exact origin of this splitting is beyond the scope of this work, and may be carried out in future studies. The transfer matrix method is utilized to extract the in-plane anisotropic complex dielectric function of several-layer ReS\textsubscript{2} in unprecedented detail, which is crucial for modelling photonic devices or further experiments that make use of its anisotropic optical properties. Our temperature dependent study reveals an anomalous temperature variation in the PL intensity which underpins the existence of low-lying dark states, indicating quasi-indirect band gap in this system. The proposed mechanism considering the dynamical processes not only explains the observed temperature variation but also estimates the phonon energies supported by previous reports. It also establishes temperature dependent study as an accessible method to probe the existence of low-lying optically dark states.

\begin{acknowledgments}
We acknowledge Dr. D. K. Goswami and his group for the AFM measurement on the sample. S. D. acknowledges SERB Ramanujan Fellowship, ISIRD grant, IIT Kharagpur and MHRD for the funding and support for this work. D. C. acknowledges CSIR for the financial assistance. We thank Dr. C. Chakraborty, and Dr. P. K. Chakraborty for their valuable comments on this work.
\end{acknowledgments}

\providecommand{\noopsort}[1]{}\providecommand{\singleletter}[1]{#1}%