\section{Introduction}
Understanding transient states in warm dense matter (WDM) is one of the grand challenges of plasma physics that is currently being tackled in a number of experimental facilities~\cite{MBRKA09, tschentscher_photon_2017, MacDonald_POP21,hoffmann_cpp_18}.
In these experiments, WDM is generated, for example, due to laser-induced shock compression~\cite{Fortov_book, graziani-book, falk_wdm}. At the foundational level, probing WDM facilitates a better understanding of astrophysical objects such as planetary interiors~\cite{AGP99,NH04, Militzer_2008, manuel, NFKR11, KKNF12} and stars~\cite{CBFS00, saumon1, Daligault_2009, Haensel}.
Furthermore, understanding WDM has great technological potential as it occurs during the fuel compression processes in inertial confinement fusion~\cite{hu_ICF}. Finally, exploring novel materials properties by driving matter through the WDM regime is an active area of research~\cite{doe-report-17, TRB08,Brongersma2015, KSM07}.

The interpretation of WDM experiments relies on a strong interplay with theory and simulation, because typical parameters like density and temperature cannot be inferred solely from the experimental data.
The use of quantum Monte Carlo (QMC) techniques has recently paved the way for determining the exact properties of interacting electrons in the uniform electrons gas under WDM conditions~\cite{dornheim_physrep_18, dornheim_jcp_19-nn, dornheim_dynamic, dornheim_prl_18, dornheim_prl_20}. However, this method is computationally expensive and does not take into account the coupling of the electrons to the ions explicitly.
Kohn-Sham density functional theory (KS-DFT)\cite{KS65} has therefore emerged as the standard method to model the electronic structure in WDM, because it includes both electron-electron and electron-ion interactions.
Due to its balance of accuracy and computational efficiency~\cite{doi:10.1146/annurev-physchem-040214-121420}, KS-DFT enables calculating materials properties such as structural and electronic transport properties under the conditions relevant to the WDM community~\cite{ISI:000177971300006, ISI:000185240600035, ISI:000227459400099, PhysRevB.77.184201, PhysRevLett.102.115701, PhysRevLett.108.091102, KDBLCSMR15}.

Here, the central quantity is the KS potential -- a mean-field potential that mimics the electron-electron interaction, in principle, exactly. The key ingredient to the KS potential is the exchange-correlation (XC) functional. In practice, the XC functional needs to be approximated. From a historical perspective, the development of XC approximations has focused on the electronic structure of molecules and solids under ambient conditions~\cite{doi:10.1063/1.4704546}. These are commonly ranked in terms of increasing accuracy and computational cost on the so-called Jacob's ladder~\cite{doi:10.1063/1.1390175}. Starting at the lowest rung with the local density approximation (LDA)~\cite{KS65}, an array of XC approximations has been developed including, for instance, the generalized gradient approximation (GGA)~\cite{PhysRevA.38.3098, PBE}, the meta-GGA~\cite{PhysRevLett.91.146401, SCAN}, and hybrid functionals~\cite{doi:10.1063/1.464913, doi:10.1063/1.472933, doi:10.1063/1.1564060}.
A key step most relevant for WDM modeling is the generalization of KS-DFT to finite temperature~\cite{mermin_65}. Based on this, several works have fleshed out the theoretical aspects of functional construction at finite temperature~\cite{PhysRevLett.107.163001, PhysRevB.84.125118, PhysRevB.93.195132, PhysRevB.93.205140}. Most recently, these have led to the construction of XC functionals that have an explicit temperature dependence \cite{PhysRevB.86.115101, PhysRevLett.112.076403, PhysRevB.88.161108, PhysRevB.88.195103, groth_prl, PhysRevB.101.245141}.

Despite these efforts, functional development has not taken into account the needs of WDM physics  to a large extent beyond the inclusion of finite-temperature effects in the electronic structure~\cite{karasiev_importance,PhysRevE.103.013210}. Some exact conditions imposed for chemistry and solid-state physics are not relevant to WDM.
For example, in contrast to solid-state physics, the surface energy has no significance in WDM, because the generated samples have no well defined surface, but rather transition to a plasma state~\cite{Olson, PhysRevE.98.063206, PhysRevE.68.036403}. Conversely, there are properties that are highly important in WDM, but have little value for solid-state physics. The most prominent example are perturbations imposed on the electronic structure with respect to the wave number $q$. An accurate description of such perturbations across a large range of $q$ is essential for WDM modeling.
Accuracy is not only required in the long wavelength regime $q<2q_F$ which is dominated by collective electronic excitations, but also at large wave numbers $q>2q_F$ where so-called single-particle effects become important~\cite{hamann_prb_20, hamann_cpp_20, zhandos_pop18, moldabekov_pop15}. Achieving a high accuracy on a wide range of wave numbers is essential, e.g., for modeling X-Ray Thomson scattering (XRTS) experiments, which is a technique of paramount importance in WDM diagnostics~\cite{GGLR2003:demonstration}.

In this paper, we therefore analyze the accuracy of common XC functionals in the WDM regime when perturbations on the electronic structure are imposed across a large range of wave numbers $q$.
To that end, we benchmark the results of KS-DFT calculations against path-integral QMC data which are considered exact within the given error bars.
Specifically, we employ four common XC functionals $-$ the LDA in the Perdew-Zunger parametrization~\cite{Perdew_LDA}, the GGAs PBE~\cite{PBE} and PBEsol~\cite{PBEsol}, and the meta-GGA SCAN~\cite{SCAN}.
We focus on the sensitivity of these XC functionals towards a perturbed electronic structure at a finite wave number $q$.
Thereby, we uncover failures of these XC approximations.
Our result strongly highlights the need for novel XC functionals that remain accurate under a perturbed electronic structure that is typically present in WDM.
Furthermore, to our knowledge, this assessment of KS-DFT against QMC data with respect to electronic density perturbations in the WDM regime has not been performed yet.
While recent developments introduced above have addressed the explicit temperature dependence of XC functionals, they do not account for the physics that stems from perturbations at finite $q$. In our analysis, we therefore separate the influence of finite temperature in the electronic states from the effect of perturbations in $q$. We follow the common practice of using the listed ground-state XC approximations and including only the implicit temperature dependence of the electronic structure in terms of a Fermi-Dirac occupation of the KS states.

The paper is organized as follows: we introduce the theoretical aspects and the simulation methods in Section~\ref{s:theory}; we present our results on benchmarking common XC approximations in Section~\ref{s:results}; we address the issue of finite size effects in Section~\ref{s:finite_size}; we assess the performance of the aforementioned XC functionals in terms of the total energy in Section~\ref{s:energy}; we provide conclusions and an outlook on future perspectives in Section~\ref{s:end}. We provide methodological details on KS-DFT and path-integral Monte Carlo (PIMC) in the Appendix.

\section{Theory and Simulation Methods}\label{s:theory}
We begin our assessment with the following Hamiltonian that imposes electronic perturbations on a warm dense uniform electron gas (UEG):
\begin{align}\label{eq:H}
\hat{H} = \hat{H}_{\rm UEG} +  \sum_{k=1}^N \sum_{i=1}^{N_p} 2\, A\cos(\hat{\bm{r}}_k \cdot \hat{\bm{q}}_i),
\end{align}
where $\hat{H}_{\rm UEG}$ denotes the standard Hamiltonian of the UEG~\cite{loos,quantum_theory,review}, $N$ the number of electrons, $N_p$ the number of external harmonic perturbations, $A$ their amplitude, and $\vec q_i$ the wave vector of the perturbations. In this work, we consider both one harmonic perturbation $N_p=1$ and the combination of two harmonic perturbations $N_p=2$. Note that we work within atomic units throughout.

It has been shown that the electronic states described by Eq.~(\ref{eq:H}) are generated in WDM experiments~\cite{PhysRevLett.125.085001}. Furthermore, Eq.~(\ref{eq:H}) has been used to rigorously examine various fundamental physical properties of the electronic structure in WDM, such as the local field correction (LFC)~\cite{PhysRevLett.75.689, PhysRevLett.69.1837, dornheim_pre17, groth_jcp17} and the non-linear response~\cite{PhysRevLett.125.085001, PhysRevLett.125.235001, dornheim2021density} that are used to describe XRTS signals. Moreover, it turns out that the UEG model provides an excellent description of many electronic properties in WDM, because the electron-ion coupling can be relatively weak~\cite{PhysRevLett.125.235001, CPP2015, GRABOWSKI2020100905}.

The Hamiltonian in Eq.~(\ref{eq:H}) is a convenient device to control the degree of inhomogeneity and the wave number of the imposed perturbations on the UEG by tuning the parameters $A$ and $\vec q$, respectively.
We analyze the performance of XC functionals with respect to electronic density perturbations by considering an amplitude range $0.02\leq A\leq 5$. This range covers regimes from weak density perturbations with $\delta n/n_0\ll1$ to strongly inhomogeneous electronic systems with $\delta n/n_0\gg 1$, where $\delta n=n-n_0$.
Additionally, we consider perturbations imposed by tuning the wave number in the range $0.843~q_F \leq q_i\leq 5.9~q_F$. Thereby, we assess the performance of XC functionals at different length scales ranging from the long wavelength regime defined by collective behavior $q<q_F$~\cite{hamann_prb_20, hamann_cpp_20} to the scale defined by the single-particle limit $q\gg q_F$~\cite{zhandos_pop18, dornheim_jcp_19-nn}. This covers the entire range of wave numbers relevant to WDM generated in experiments.

WDM conditions prevail when we consider matter at solid density and a temperature $T \sim T_F$ \cite{bonitz_pop_08, moldabekov_pre_18}, where $T_F=E_F/k_B$ denotes the Fermi temperature defined in terms of the Fermi energy $E_F$ and the Boltzmann constant $k_B$. In our analysis, we consider the densities $r_s=2$ and $r_s=6$ and a temperature $T=T_F$, where $r_s=a/a_B$ defines the number density of electrons which is given as the ratio between the Wigner-Seitz radius $a$ and the first Bohr radius $a_B$. At $T\lesssim T_F$, the parameter $r_s$ also characterizes electronic non-ideality~\cite{ott_epjd18, Moldabekov_cpp_21}. Therefore, it is used also as a coupling parameter.
These parameters are typically encountered in experiments with laser-driven and  shock-compressed solid targets~\cite{LANDEN2001465, PhysRevLett.98.065002, PhysRevLett.109.065002}.

The KS-DFT calculations were performed with the GPAW code~\cite{GPAW1, GPAW2, ase-paper, ase-paper2}, which is a real-space implementation of the projector augmented-wave method. Details of the KS-DFT simulation parameters such as the number of $\vec k$-points and number of bands are given in the Appendix.

We provide unassailable data for benchmarking our KS-DFT results by carrying out direct path-integral QMC calculations based on Eq.~(\ref{eq:H}) without any restrictions on the nodal structure of the thermal density matrix. Therefore, the calculations are computationally expensive due to the fermion sign problem~\cite{dornheim_sign_problem, dornheim2021fermion}, but \emph{exact} within the given Monte Carlo error bars. The corresponding simulation details, such as the number of imaginary-time propagators, can be found in the Appendix~\ref{sec:app2}.

\begin{figure}
\center
\includegraphics{density_plot_A01.jpg}
\caption{ \label{fig:den_rs2} The electronic density distribution along the perturbation direction for two different amplitudes $A$ and increasing wave number $q$ at $r_s=2$ and $T=T_F$.
}
\end{figure}
\begin{figure*}
\center
\includegraphics{dnEr_fig_1qA0_02.jpg}
\caption{\label{fig:den_err1}
Relative deviation in the density $\Delta n/{\rm max}\{\delta n\} 100~\%$ between the KS-DFT data and the reference QMC data at $r_s=2$. Column-wise from left to right: $A=0.02,~0.1,~0.5$ and $A=1$.
Row-wise from top to bottom: $q_1=q_{\rm min},~2 q_{\rm min}$, and $q_1=3 q_{\rm min}$.}
\end{figure*}

\section{Results}\label{s:results}

\subsection{Single Harmonic Perturbation}
Let us begin our analysis by considering the case of a single harmonic perturbation, i.e., $N_p=1$ in Eq.~(\ref{eq:H}).
First, we investigate the electronic density distribution at a density (coupling) parameter $r_s=2$ for a range of perturbation amplitudes and wave numbers. Then, we consider a stronger coupling regime at $r_s=6$. These parameters are of particular interest in WDM experiments.

\subsubsection{Metallic Density, $r_s=2$}
The density distribution along the direction of the density perturbation is shown for two amplitudes $A=0.1$ and $A=1$ and increasing wave number in Fig.~\ref{fig:den_rs2}.
The top panel displays the results for $q_1=q_{\rm min}=0.843~q_F $, the middle panel for $q_1=2 q_{\rm min}$, and the bottom panel for $q_1=3 q_{\rm min}$. On the scale of the total density, all considered XC functionals are in overall agreement with each other and with the QMC data.

Let us now consider the actual deviation of the KS-DFT data from the reference QMC data in closer detail.
To that end, we consider the relative density deviation $\Delta n/{\rm max}\{\delta n\}$ between the KS-DFT data and the reference QMC data, where ${\rm max}\{\delta n\}$ is the maximum deviation of the QMC data from the mean density. We use ${\rm max}\{\delta n\}$ for the analysis of the KS-DFT results because the physically important quantity in the case of a weak perturbation $\delta n/n_0\ll1$ is the deviation of the density from $n_0$ rather than the total density $n$ itself. Indeed, $\delta n$ defines the density response of the system which is a cornerstone of linear response theory describing all related physical properties of electrons in equilibrium.
More specifically, when the external field defined via Eq.~(\ref{eq:H}) is weak, the density perturbation $\delta n=2A \chi(\vec q) \cos(\vec q\cdot \vec r)$ defines the static response function $\chi(\vec q)$. Therefore, $\Delta n/{\rm max}\{\delta n\}=\Delta \chi/\chi_{\rm QMC}$ measures the error in the response function.
On the other hand, to keep the present analysis general, we provide the values of ${\rm max}\{\delta n\}$ along with $\Delta n/{\rm max}\{\delta n\}$. This allows a simple conversion of data from $\Delta n/n_0$ or $\Delta n/n$.

In Fig.~\ref{fig:den_err1}, we show the relative difference between the KS-DFT data and the reference QMC data at $r_s=2$ for the amplitudes $A=0.02,~0.1,~0.5,~1$ (from left to right) and for the wave numbers $q_1=q_{\rm min},~2 q_{\rm min},~3 q_{\rm min}$ (from top to bottom). Additionally, the corresponding largest absolute value of the discrepancy is given in Table \ref{Table1}.

\begin{table*}[]
\caption{The performance of common XC functionals in terms of the relative density deviation $\Delta n/{\rm max}\{\delta n\} 100~\%$. A single harmonic perturbation at a fixed density $r_s=2$ (metallic density) and varying perturbation amplitude $0.02\leq A\leq 1$ and wave number $q_{\rm min}\leq q\leq 3 q_{\rm min}$ is considered, where $q_{\rm min}=0.843~q_F$. The largest absolute values of the deviation are listed in this table.}
\vspace{0.5cm}
\label{Table1}
\begin{tabular}{ c m{1em} ccc m{1em} ccc m{1em} ccc m{1em} ccc }
\hline
\hline
\multirow{2}{*}{} & & \multicolumn{3}{c}{$\bf A=0.02$} & & \multicolumn{3}{c}{$\bf A=0.1$} & & \multicolumn{3}{c}{$\bf A=0.5$} & & \multicolumn{3}{c}{$\bf A=1.0$} \\
\hline
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf  3 q_{\rm min}$
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf 3 q_{\rm min}$
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf 3 q_{\rm min}$
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$  & $\bf 3 q_{\rm min}$\\
 \hline
{\bf LDA}    & & 6.59 & 2.01 & 4.07 & & 2.20 & 1.55 & 4.86 & & 2.37 & 1.57 & 4.22 & & 2.49 & 2.10 & 1.42\\
\hline
{\bf PBE}    & & 6.59 & 2.02 & 4.09 & & 2.20 & 1.72 & 5.34 & & 2.41 & 3.20 & 7.50 & & 2.52 & 3.36 & 2.18\\
\hline
{\bf PBEsol} & & 6.72 & 2.08 & 2.73 & & 2.40 & 1.40 & 5.60 & & 2.56 & 2.28 & 5.52 & & 2.72 & 2.71 & 1.92\\
\hline
{\bf SCAN}   & & 5.66 & 5.16 & 7.09 & & 0.95 & 4.24 & 8.75 & & 0.69 & 4.33 & 6.27 & & 0.86 & 2.29 & 1.49\\
\hline
\hline
\end{tabular}
\end{table*}

\begin{figure}
\center
\includegraphics{7qA1.jpg}
\caption{\label{fig:den_err2}
Same as in Fig.~\ref{fig:den_err1}, but for $q_1=~7 q_{\rm min}$ at $A=1$ (top ) and $A=5$ (bottom).}
\end{figure}

\begin{figure}
\center
\includegraphics{density_plot_A5.jpg}
\caption{\label{fig:density_A5}
The electronic density distribution along the perturbation direction for $A=5$ and $q_1=~7 q_{\rm min}$ at $r_s=2$.}
\end{figure}

Our assessment of the results presented in Fig.~\ref{fig:den_err1} and Table~\ref{Table1} can be summarized as follows.
At $A=0.02$ and $q_1=q_{\rm min}$, the results computed with the LDA, PBE, and PBEsol functionals have about the same level of accuracy. They are consistent with the QMC reference densities to within $7~\%$. With a maximum deviation of $5.66~\%$ from the QMC data, the SCAN functional yields slightly more accurate results. Note that here the statistical uncertainty of the QMC results is of the same order as the relative difference between the KS-DFT and QMC data (see the top left corner in Fig.~\ref{fig:den_err1}). Upon increasing the perturbation strength up to $A=1.0$, the relative error of the KS-DFT obtained from the LDA and GGA functionals remains less than $\leq 3~\%$, while the use of SCAN provides a remarkable accuracy better than $1~\%$. Therefore, the KS-DFT calculations using SCAN are virtually exact at $q_1=0.843~q_F $ for both weak and strong perturbations.
Contrarily, the accuracy provided by the SCAN functional is not maintained with an increase of the perturbation in terms of its wave number to $q_1=2q_{\rm min}$ and $q_1=3 q_{\rm min}$. Across the perturbation amplitdes $0.02\leq A\leq 0.5$ ($\delta n/n_0\leq 0.244$), the relative differences for calculations using SCAN are about $4~\%$ and $8~\%$. The LDA and PBEsol functionals provide an accuracy of about $2~\%$ ($5~\%$) at $q_1=2 q_{\rm min}$ ($q_1=3 q_{\rm min}$). The PBE functional provides a comparable accuracy to LDA and PBEsol at $q_1=2 q_{\rm min}$, but becomes less accurate at $q_1=3 q_{\rm min}$ and $A=0.5$ reaching an error of $7.5~\%$.
In the strong-perturbation regime ($A=1.0$ and ${\rm max}\{\delta n\}>1.3~n_0$), all functionals LDA, PBE, PBEsol functionals now provide about the same level of accuracy when compared to the QMC data.
\textit{In summary, overall the LDA and PBEsol show a more robust performance compared to SCAN and PBE for all perturbation amplitudes and relatively small wave numbers at the typical mass density of metals.}

Next in Fig.~\ref{fig:den_err2}, we assess the performance of the XC functionals in the limit of large perturbation wave numbers by setting $q_1=7 q_{\rm min}=5.9~q_F$, but keeping the mass density metallic ($r_s=2$). We consider the perturbation amplitudes $A=1$ (with ${\rm max}\{\delta n\}= 0.28~n_0$) and $A=5$ (with ${\rm max}\{\delta n\}= 1.268~n_0$).
At $A=1$, we find that for LDA, PBEsol, and SCAN, the relative difference between the KS-DFT data and the QMC data is less than $2.5~\%$, while the PBE functional yields an error of $5~\%$.
When we increase the perturbation amplitude to $A=5$, the KS-DFT results from all considered functionals are in agreement with each other, but exhibit a strong disagreement with the QMC data of up to about $12~\%$. In this extreme regime, the deviation of the KS-DFT data from the QMC results is significant even on the scale of the total density $n$. This is illustrated in Fig.~\ref{fig:density_A5}.
\textit{In summary, all tested XC functionals fail dramatically in yielding an accurate electronic density for $q\gg q_F$ and ${\rm max}\{\delta n\}/n_0>1$ at metallic density.  We argue that a large number of other XC functionals that are derived from LDA, GGA, and meta-GGA classes are afflicted with the same limitation.}

\begin{figure}
\center
\includegraphics{density_plot2.jpg}
\caption{\label{fig:density_rs6}
The electronic density distribution along the perturbation direction for two different amplitudes ($A=0.05$ and $A=0.1$) at $r_s=6$.}
\end{figure}
\begin{figure*}
\center
\includegraphics{rs6dnEr_fig_1qA0_02.jpg}
\caption{\label{fig:den_err_rs6}
Relative deviation in the density $\Delta n/{\rm max}\{\delta n\} 100~\%$ between the KS-DFT data and the reference QMC data at $r_s=6$. Column-wise from left to right: $A=0.01,~0.05$ and $A=0.1$.
Row-wise from top to bottom: $q_1=q_{\rm min},~2 q_{\rm min}$, and $q_1=3 q_{\rm min}$.}
\end{figure*}

\subsubsection{Strong Coupling, $r_s=6$}
Next, we investigate strongly correlated electronic systems where $r_s=6$. Such low densities can be realized, for example, in evaporation experiments~\cite{Zastrau}.
From a theoretical perspective, these conditions are particularly challenging due to the substantial impact of electronic XC effects~\cite{low_density1,low_density2} on physical observables like the electrical conductivity or the density.

Again, we consider weak as well as strong perturbation amplitudes corresponding to a density enhancement in the range $0.034~n_0\leq \delta n \leq 2.044~n_0$.
The electronic density distribution along the direction of the perturbation is shown in Fig.~\ref{fig:density_rs6} for $q_1=q_{\rm min},~2 q_{\rm min}$, and $q_1=3 q_{\rm min}$ (from top to bottom) at $A=0.05$ (left) and $A=1.0$ (right).
The agreement between the KS-DFT data and QMC data is excellent when $q_1=q_{\rm min}$. With increasing wave number, the accuracy decreases. For example, the errors are larger than in the previous case where a metallic density is considered (see Fig.~\ref{fig:den_rs2}).

We further delineate the relative difference between the KS-DFT data and the reference QMC results in Fig.~\ref{fig:den_err_rs6}.
There, the relative difference is illustrated for the amplitudes $A=0.01,~0.05$ and $A=0.1$ (from left to right) and the wave number $q_1=q_{\rm min},~2 q_{\rm min}$, and $q_1=3 q_{\rm min}$ (top to bottom) at $r_s=6$. The largest absolute values of the differences are listed in Table~\ref{Table2}.

The main conclusions of our assessment in this parameters range are as follows.
At $q_1=q_{\rm min}$ and all perturbation amplitudes, the SCAN functional yields the most reliable results with an accuracy better than $1.64~\%$. The other tested  functionals provide an accuracy better than $3.13~\%$.
In contrast to that, when we increase the wave number to $q_1=2 q_{\rm min}$, SCAN performs much worse with a maximum deviation of $12.58~\%$ at $A=0.01$ and $9.97~\%$ at $A=0.05$. Both the LDA and PBEsol functional yield an accuracy better than $3.5~\%$. Similarly, PBE yields an error of $4.02~\%$ and $5.08~\%$ at $A=0.05$ and $A=0.1$, respectively.
However, irrespective of the wave number, the SCAN functional provides the most accurate results with the maximum deviation of $2.63~\%$ at a stronger perturbation amplitude $A=0.1$.
Finally, increasing the wave number of the perturbation to $q_1=3 q_{\rm min}$ renders KS-DFT data less accurate with a maximal deviation in the range of about $8~\%$ and $12~\%$ (for $0.01\leq A\leq 0.1$) when LDA and PBE are used. The PBEsol functional is in better agreement with the QMC results showing a maximum deviation in the range between $5.63~\%$ and $9.45~\%$. The SCAN functional provides very low accuracy when the perturbation is weak($A=0.01$), with a deviation of almost $30~\%$. This improves in the case of strong perturbation $A=0.1$, where the maximum deviation is $6.2~\%$.
\textit{In summary, the overall performance of the considered XC functionals is worse when strongly coupled electronic systems are considered. We observe failures of the functionals, in particular, when the wave number of the perturbation increases.}

\begin{table*}[]
\caption{The performance of common XC functionals in terms of the relative density deviation $\Delta n/{\rm max}\{\delta n\} 100~\%$. A single harmonic perturbation at a fixed density $r_s=6$ (strong coupling) and varying perturbation amplitude $0.01\leq A\leq 0.1$ and wave number $q_{\rm min}\leq q\leq 3 q_{\rm min}$ is considered, where $q_{\rm min}=0.843~q_F$. The largest absolute values of the deviation are listed in this table.}
\vspace{0.5cm}
\label{Table2}
\begin{tabular}{ c m{1em} ccc m{1em} ccc m{1em} ccc}
\hline
\hline
\multirow{2}{*}{} & & \multicolumn{3}{c}{$\mathbf{A=0.01}$} & & \multicolumn{3}{c}{$\mathbf{A=0.05}$} & & \multicolumn{3}{c}{$\mathbf{A=0.1}$} \\
\hline
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf  3 q_{\rm min}$
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf 3 q_{\rm min}$
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf 3 q_{\rm min}$\\
 \hline
{\bf LDA}    & & 2.37 &  1.18 & 11.81 & & 2.66 & 2.29 & 10.60 & & 3.01 & 3.44 &  7.85\\
\hline
{\bf PBE}    & & 2.37 &  1.18 & 11.84 & & 2.61 & 4.02 & 15.57 & & 3.04 & 5.08 & 10.68\\
\hline
{\bf PBEsol} & & 2.63 &  2.04 &  5.63 & & 2.80 & 1.16 &  9.45 & & 3.13 & 3.48 &  8.86\\
\hline
{\bf SCAN}   & & 0.91 & 12.58 & 29.85 & & 1.55 & 9.97 & 14.02 & & 1.64 & 2.63 &  6.20\\
\hline
\hline
\end{tabular}
\end{table*}

\subsection{Double Harmonic Perturbation}
Furthermore, we assess the accuracy of the considered XC functionals when more complex perturbations are applied. To that end, we consider a double harmonic perturbation $N_p=2$ with $q_1=q_{\rm min}$ and $q_2=2 q_{\rm min}$. This allows us to check if the observed poor performance of the XC functionals manifests itself when the perturbations of different wave numbers are superimposed. Again, we consider both metallic densities ($r_s=2$) and strongly coupled systems ($r_s=6$).
When $r_s=2$, we set the perturbation amplitude to $A=0.1$ resulting in a density perturbation of the order of $0.1~n_0$. When $r_s=6$, we set $A=0.01$ leading to a similar density perturbation of the order of $0.1~n_0$.

The resulting electronic density distributions are shown in Fig.~\ref{fig:density_two_A}. At $r_s=2$, the KS-DFT results show good agreement with the QMC data when gauged on the scale of the total density. At $r_s=6$, the SCAN functional deviates from the QMC data significantly, while the results obtained from the other XC functionals are indistinguishable from the QMC data on this scale.
We take a closer look at the performance of the various XC functionals in Fig.~\ref{fig:density_err_two_A}, where we use the relative deviation in the density $\Delta n/{\rm max}\{\delta n\}$ between the KS-DFT data and the reference QMC results. The corresponding largest absolute values of the differences are given in Table~\ref{Table3}.

The results in Fig.~\ref{fig:density_err_two_A} and Table \ref{Table3} show that LDA, PBE, and PBEsol provide an accuracy better than $2.3~\%$ for both $r_s=2$ and $r_s=6$. Contrarily, the densities computed using SCAN have a maximum deviation of $4.31~\%$ and $12.05~\%$ at  $r_s=2$ and $r_s=6$, respectively. These numbers are similar to the deviations observed in the case of a single harmonic perturbation at the wave number $2 q_{\rm min}$.
This provides a strong indication about the general applicability of the present findings to other systems, as any external potential can be expressed as a superposition of harmonic perturbations in reciprocal space.

\begin{table}[]
\caption{The largest absolute value of $\Delta n/{\rm max}\{\delta n\} 100~\%$ for $r_s=2$ with $A=0.1$  and for $r_s=6$ with $A=0.01$. The wave numbers of the double harmonic perturbation in Eq.~(\ref{eq:H}) are $q_1=q_{\rm min}$ and $q_2=2 q_{\rm min}$. }
\vspace{0.5cm}
\label{Table3}
\begin{tabular}{ c c c }
\hline
\hline
$\mathbf{r_s}$ & {\bf 2.0} & {\bf 6.0} \\
\hline
{\bf LDA}      &  1.89   &  1.66   \\
\hline
{\bf PBE}      &  2.26   &  1.66   \\
\hline
{\bf PBEsol}   &  1.82   &  2.27   \\
\hline
{\bf SCAN}     &  4.31   &  12.05  \\
\hline
\hline
\end{tabular}
\end{table}

\begin{figure}
\center
\includegraphics{density_q1q2_A0_01.jpg}
\caption{\label{fig:density_two_A}
The density distribution along the perturbation direction for $r_s=2$ with $A=0.1$ (left) and for $r_s=6$ with $A=0.01$ (right). The wave numbers of the double harmonic perturbation in Eq.~(\ref{eq:H}) are $q_1=q_{\rm min}$ and $q_2=2 q_{\rm min}$. }
\end{figure}
\begin{figure}
\center
\includegraphics{dnerr_q1q2_A0_01.jpg}
\caption{\label{fig:density_err_two_A}
Relative deviation in the density $\Delta n/{\rm max }~\delta n \times 100~\%$ between the KS-DFT data and the reference QMC data for $r_s=2$ with $A=0.1$ (left) and for $r_s=6$ with $A=0.01$ (right). The wave numbers of the double harmonic perturbation in Eq.~(\ref{eq:H}) are $q_1=q_{\rm min}$ and $q_2=2 q_{\rm min}$.}
\end{figure}

\section{Finite size effects}\label{s:finite_size}

\begin{figure}
\center
\includegraphics{N20.jpg}
\caption{\label{fig:den_N20}
The electronic density distribution along the perturbation direction for two different values of wave numbers $q=q_{\rm min}$ (top) and $q=3q_{\rm min}$ (bottom) at $A=0.1$ and $r_s=2$.}
\end{figure}

\begin{figure}
\center
\includegraphics{N20dif.jpg}
\caption{\label{fig:dif_N20}
Relative deviation in the density $\Delta n/{\rm max}\{\delta n\} 100~\%$ between the KS-DFT data and the reference QMC data for  $q=q_{\rm min}$ (top) and $q=3q_{\rm min}$ (bottom) at $r_s=2$ and $A=0.1$.}
\end{figure}

\begin{figure}
\center
\includegraphics{N34.jpg}
\caption{\label{fig:den_N34}
The electronic density distribution along the perturbation direction for two different values of wave numbers $q=q_{\rm min}$ (top) and $q=3q_{\rm min}$ (bottom) at $A=0.1$ and $r_s=2$.}
\end{figure}

\begin{figure}
\center
\includegraphics{N34dif.jpg}
\caption{\label{fig:dif_N34}
Relative deviation in the density $\Delta n/{\rm max}\{\delta n\} 100~\%$ between the KS-DFT data and the reference QMC data for  $q=q_{\rm min}$ (top) and $q=3q_{\rm min}$ (bottom) at $r_s=2$ and $A=0.1$.}
\end{figure}

The results presented so far were computed for $14$ electrons in a simulation cell. Naturally, it is important to inquire whether the large deviations of the KS-DFT results from the QMC data observed at $2q_{\rm min}$ and $3q_{\rm min}$ compared to the case with $q_1=q_{\rm min}$ are indeed genuine results or an artifact due to finite size effects.
To exclude this possibility, we consider $N=20$ and $N=34$ electrons in the simulation cell at $r_s=2$ and $A=0.1$.
The density distribution along the direction of the density perturbation for the case with $20$ electrons is shown in Fig.~\ref{fig:den_N20}. The top panel displays the results for $q_1=q_{\rm min}=0.843~q_F $ and the bottom panel for $q_1=3 q_{\rm min}$.
We observe that on the scale of the total density, the KS-DFT results are in agreement with the QMC data.
The corresponding relative difference between the KS-DFT results and the reference QMC data is displayed in Fig.~\ref{fig:dif_N20}, where the top panel is for $q_1=q_{\rm min}$ and the bottom panel is for $3 q_{\rm min}$. Fig.~\ref{fig:dif_N20} shows that at $q_1=q_{\rm min}$ the relative difference between KS-DFT results obtained using the SCAN functional and the QMC data is within the given uncertainty range, whereas it increases to up to $8\%$ with an increase in the wave number of the perturbation to $3q_{\rm min}$. The other considered XC functionals lead to slightly worse agreement at $q_1=q_{\rm min}$, with about $2\%$ maximum relative difference, and to much better results at $q_1=3q_{\rm min}$, with a maximum relative difference of about $4\%$.

Next, we further increase the number of electrons in the simulation cell to $34$. The density distribution for $34$ particles along the direction of the density perturbation is shown in Fig.~\ref{fig:den_N34} at  $q_1=q_{\rm min}$ (top panel) and at $q_1=3q_{\rm min}$ (bottom panel).
Again, on the scale of the total density we observe agreement of KS-DFT results with the QMC data. The corresponding relative differences, shown in Fig.~\ref{fig:dif_N34}, confirm this assessment. At $q_1=q_{\rm min}$ (top panel) the KS-DFT results agree with the QMC data within the given numerical uncertainty. At $q_1=3q_{\rm min}$ (the bottom panel in Fig.~\ref{fig:dif_N34}), KS-DFT results obtained with the SCAN functional exhibit a relative deviation from QMC data of about  $8\%$. The other XC functionals yield a maximum relative difference of about $4\%$. Note the statistical uncertainty of the QMC data at $N=34$ is noticeably larger due to the infamous fermion sign problem. See Appendix for more details and Ref.~\onlinecite{dornheim_sign_problem} for a topical review of this issue.

This analysis of the densities and relative differences between KS-DFT results and QMC data for $N=20$ and $N=34$ electrons clearly shows that our conclusions at $N=14$ electrons are not affected by finite size effects.

\section{Energy perturbation}\label{s:energy}

\begin{figure}
\center
\includegraphics{JCP_fig_energy.jpg}
\caption{\label{fig:energy}
Relative deviation in the total energy from UEG ($A=0$) due to external perturbation with amplitude $A$ (see Hamiltonian (\ref{eq:H})) at $r_s=2$. From top to bottom: $q_1=q_{\rm min},~2 q_{\rm min}$, and $q_1=3 q_{\rm min}$.}
\end{figure}

\begin{figure}
\center
\includegraphics{JCP_fig_energyrs6.jpg}
\caption{\label{fig:energy_rs6}
Relative deviation in the total energy from UEG ($A=0$) due to external perturbation with amplitude $A$ (see Hamiltonian (\ref{eq:H})) at $r_s=6$.
Row-wise from top to bottom: $q_1=q_{\rm min},~2 q_{\rm min}$, and $q_1=3 q_{\rm min}$.}
\end{figure}

Finally, we analyze the correlation of the revealed inaccuracies in the densities with the total energy.
To that end, we focus on the change in the total energy compared to the unperturbed case. This is quantified by the reduced perturbation in the total energy
\begin{equation}\label{eq:dE}
    \delta \widetilde{E}= \frac{\left|E(A)-E(A=0)\right|}{\left|E_{\rm QMC}(A=0)\right|},
\end{equation}
where $E(A)$ is the total energy of electrons in the presence of the external perturbation with the amplitude $A$, $E(A=0)$ is the total energy of the unperturbed system, and $E^{\rm QMC}(A=0)$ is the total energy of the unperturbed system from QMC simulations.

The results for $\delta \widetilde{E}$ at $r_s=2$ and at  $r_s=6$ for different amplitudes and wave numbers are presented in Fig.~\ref{fig:energy} and Fig.~\ref{fig:energy_rs6}, respectively. In these figures we also plot $\delta \widetilde{E}\sim A^2$ using the data point at $A=0.1$ ($A=0.01$) for $r_s=2$ ($r_s=6$). This dependence is known to be exact in the limit of small $A$ in the ground state at $T\ll T_F$ and is referred to as the stiffness theorem \cite{giuliani2005quantum}. At finite temperature, the stiffness theorem applies to the total free energy. Figs.~\ref{fig:energy} and \ref{fig:energy_rs6} can be understood as an empirical indication of validity of stiffens theorem for the total energy at finite temperature.
The deviation of the simulation data from the $A^2$ dependence is related to the emergence of the higher-order non-linear effects.

In Fig.~\ref{fig:energy}, in addition to KS-DFT results within LDA, PBE, PBEsol, and SCAN, we also show results computed using null XC functional, i.e., the electron-electron interaction approximated only using the mean Hartree field. First of all, in Fig.~\ref{fig:energy} we clearly see that neglecting the XC energy leads to less accurate results compared to the QMC data at amplitudes $A=0.1$, $A=0.5$ and $A=1.0$. At $A=0.02$, the change in the energy due to the external field is comparable with statistical uncertainty. This prevents us from quantifying the accuracy of the KS-DFT results for a reduced perturbation in the energy at that perturbation strength.
At $A=0.1$, all KS-DFT data points, except the null XC functional case, are in agreement with the QMC data within the given statistical uncertainty.
Furthermore, in Table \ref{Table_energy}, we provide the relative difference of the reduced perturbation in the total energy, $\Delta \widetilde E$, which is defined by

\begin{equation}\label{eq:DE}
    \Delta \widetilde E ~[\%]=\frac{|\delta \widetilde{E}_{\rm QMC}-\delta \widetilde{E}_{\rm DFT}|}{\delta \widetilde{E}_{\rm QMC}} \times 100\ .
\end{equation}

\begin{table}[]
\caption{The performance of common XC functionals in terms of the relative deviation in the total energy as defined by Eq.~(\ref{eq:DE}). A single harmonic perturbation at a fixed density $r_s=2$ and  perturbation amplitudes $A=0.5$ and $A=1.0$ and wave numbers $q_{\rm min}\leq q\leq 3 q_{\rm min}$ is considered, where $q_{\rm min}=0.843~q_F$. The presented data has about $1\%$ uncertainty.}
\vspace{0.5cm}
\label{Table_energy}
\begin{tabular}{ c m{1em} ccc m{1em} ccc}
\hline
\hline
\multirow{2}{*}{} & & \multicolumn{3}{c}{$\mathbf{A=0.5}$} & & \multicolumn{3}{c}{$\mathbf{A=1.0}$} \\
\hline
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf  3 q_{\rm min}$
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf 3 q_{\rm min}$\\
 \hline
{\bf LDA}    & & 3.05  &  3.62  & 2.52 & & 3.63 & 0.70 & 2.29 \\
\hline
{\bf PBE}    & & 2.93&  5.84  & 0.75 & & 3.31 & 0.47 & 4.24 \\
\hline
{\bf PBEsol} & & 3.31 &  4.24 &  1.77 & &  3.62 & 1.36 &  3.23 \\
\hline
{\bf SCAN}   & & 0.54 & 9.92 & 5.56 & & 1.11 & 6.61 & 0.52 \\
\hline
{\bf Null XC}   & & 15.44 & 8.8 & 11.43 & & 12.84 & 7.47 & 2.98 \\
\hline
\hline
\end{tabular}
\end{table}

\begin{table}[]
\caption{The performance of common XC functionals in terms of the relative deviation in the total energy as defined by Eq.~(\ref{eq:DE}). A single harmonic perturbation at a fixed density $r_s=6$ and  perturbation amplitudes $A=0.05$ and $A=0.1$ and wave numbers $q_{\rm min}\leq q\leq 3 q_{\rm min}$ is considered, where $q_{\rm min}=0.843~q_F$. The presented data has about $1\%$ uncertainty. }
\vspace{0.5cm}
\label{Table_energy_rs6}
\begin{tabular}{ c m{1em} ccc m{1em} ccc}
\hline
\hline
\multirow{2}{*}{} & & \multicolumn{3}{c}{$\mathbf{A=0.05}$} & & \multicolumn{3}{c}{$\mathbf{A=0.1}$} \\
\hline
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf  3 q_{\rm min}$
& & $\bf q_{\rm min}$   & $\bf 2 q_{\rm min}$ & $\bf 3 q_{\rm min}$\\
 \hline
{\bf LDA}    & &  1.16 &  1.95 & 6.91& & 1.96 & 0.26 & 5.25\\
\hline
{\bf PBE}    & & 1.16 &  0.1  & 10.84& & 1.35 & 2.19 & 8.81 \\
\hline
{\bf PBEsol} & & 1.57 &  3.88 &  4.14 & &  1.80 & 0.55 &  5.80 \\
\hline
{\bf SCAN}   & & 0.43 & 17.80& 30.93 & & 3.07 & 3.15 & 7.42 \\
\hline
\hline
\end{tabular}
\end{table}

From Table \ref{Table_energy} we infer that in the case of null XC functional, $\Delta \widetilde E$ is about $15\%$ ($9\%$) and $13\%$ ($7\%$) for $A=0.5$ and $A=1.0$ at $q=q_{\rm min}$ ($q=2q_{\rm min}$). This large difference mainly stems from an inaccurate density due to neglecting the XC functional (see Appendix). This clearly demonstrates the importance of including the XC functional in the considered KS-DFT calculations.
Also, as shown in Table \ref{Table_energy}, the LDA, PBE, and PBEsol functionals show similar performance.  The fact that SCAN is more accurate for $\Delta \widetilde E$ at $q_1=q_{\rm min}$ and less accurate at $q_1=2q_{\rm min}$ correlates with the previous assessment of the densities.

In the case of $r_s=6$ (Fig.~\ref{fig:energy_rs6}), the dependence of $\delta \widetilde{E}$ on $A$ is similar to that at $r_s=2$. At $A=0.01$, the KS-DFT results computed using different XC functionals is in agreement with the QMC data within the given statistical uncertainty. At $A=0.05$ and $A=0.1$, the SCAN results tend to overestimate $\delta \widetilde{E}$ with increasing wave number from $q_1=q_{\rm min}$ to $q_1=2q_{\rm min}$ and $q_1=3q_{\rm min}$. The LDA, PBE, and PBEsol results are very similar across all considered $A$ and $q_1$ values.

In Table~\ref{Table_energy_rs6} we list $\Delta \widetilde E$ as defined by Eq.~(\ref{eq:DE}) in order to better delineate the differences. Similar to $r_s=2$, also at $r_s=6$ the SCAN results are most accurate at $q_1=q_{\rm min}$. Contrarily, at $q_1=2q_{\rm min}$ and $q_1=3q_{\rm min}$, the largest deviation from the QMC data is exhibited by the SCAN results. This reflects the fact that there are more inaccuracies in the density when SCAN is used over the other XC functionals. The KS-DFT results for $\delta \widetilde{E}$ computed using LDA, PBE, and PBEsol have similar values at $q_1=q_{\rm min}$ and $q_1=2q_{\rm min}$. At $q_1=3q_{\rm min}$, the PBE results are less accurate by a few percent than the LDA and PBEsol results. Overall, by comparing data in Table \ref{Table2} and in Table \ref{Table_energy_rs6}, we observe that the variations in $\Delta \widetilde E$ reflect those in $\Delta n/{\rm max}\{\delta n\}$.

\section{Conclusions and Outlook}\label{s:end}

We benchmarked the performance of KS-DFT based on the LDA, PBE, PBEsol, and SCAN XC functionals against exact QMC data in the WDM regime. Our assessment revealed a set of conditions for the successful simulation of WDM with KS-DFT at QMC level accuracy. Our comparative analysis unambiguously demonstrates when the considered XC functionals fail to correctly describe the electronic density.

We found that the KS-DFT results are sufficiently accurate for small wave numbers of the density perturbation, $q<q_F$. In particular, using the SCAN functional yields an excellent agreement with the QMC reference data. However, with increasing wave number $q>q_F$ (as tested for $q=1.686q_F$ and $q=2.529q_F$), the SCAN functional performs much worse than LDA, PBE, and PBEsol. This is somewhat surprising, because the SCAN functional is ostensibly designed to be superior over LDA and GGA functionals.
In contrast to that, LDA and PBEsol show a robust performance with an accuracy better than $4~\%$ at $q\leq 1.686~q_F$ for both $r_s=2$ and $r_s=6$.

As a key finding of our assessment, we highlight that the overall performance of the considered XC functionals deteriorates upon increasing the wave number of the density perturbation.
At $r_s=2$ and $q= 2.529~q_F$, LDA and PBEsol still provide an accuracy with an error of less than $6~\%$ in $\delta n=n-n_0$.
At $r_s=6$, the same is valid for PBEsol in the regime of weak perturbations, $\delta n/n_0\ll 1$. Other XC functionals essentially fail at $r_s=6$, $\delta n/n_0< 1$, and $q= 2.529~q_F$. At the largest considered wave number, $q= 5.9~q_F$, all considered XC functionals yield errors of less than $6~\%$ if $\delta n/n_0\ll 1$ for metallic densities ($r_s=2.0)$.
Finally, we highlight the failure of the considered XC functionals in the regime of strong perturbations, $\delta n/n_0>1$, and large wave-numbers, $q= 5.9~q_F$, where they exhibit a maximum deviation of about $10~\%$. We also showed that the reported errors in the density correlate with the errors in the total energy.

Based on the performed analysis we can formulate the following general recommendations for using XC functionals within the typical WDM regime at temperatures $T\leq T_F$ and when the total density is of interest:  When characteristic wave numbers $q<q_F$ are of interest, the SCAN functional is the most reliable choice, it provides accuracy at the level of QMC. For a wider range of wave numbers $q\leq 5.9~q_F$ (at $r_s=2$) and $q\leq 1.686~q_F$ (at $r_s=6$), the LDA and PBEsol functionals should be used if $\delta n/n_0<1$, because they provide consistent results with a relative error not exceeding a few percent.

One can understand the performance of the considered XC functionals by recalling that, in the case of a weak perturbation, the density deviation from the mean value is defined by the static response function $\chi(q)$ since $\delta n(r)=2A\chi(\vec q)\cos{\vec q\cdot\vec r}$. Within linear response theory of the UEG, XC effects are included in terms of the local field correction $G(\vec q)$ via the relation $\chi(\vec q)^{-1}=\chi_0(\vec q)^{-1}-\frac{4\pi}{q^2}\left(1- G(\vec q)\right)$, where $\chi_0(\vec q)$ is the response function of ideal electron gas, i.e., the non-interacting UEG. In this case, the XC functional is $K_{_\mathrm{{\scriptscriptstyle xc}}}(\vec q)=-\frac{4\pi}{q^2} G(\vec q)$. In the UEG, the LDA functional reduces to the long wavelength limit $q\ll 2q_F$ of $K_{_\mathrm{{\scriptscriptstyle xc}}}$. Clearly, this behavior can be attributed to the failure of the LDA at perturbation wave numbers $q>2q_F$. The same is valid for PBE and PBEsol, where gradient corrections to exchange and correlation exactly cancel each other in the limit of the UEG.
On the other hand, when the density perturbation is strong, the characteristic wave number of the density perturbation defined by the local gradient of the density is larger than the wave number of the external harmonic perturbation. In this regime, regions with a strong density localization are formed. These lead to an increase of the kinetic energy and reduce the importance the wave-number dependence in XC functionals. This may explain the better performance of the considered XC functionals when the perturbation amplitude is increased.
Additionally, we note that the SCAN functional is designed to yield PBE-like results in the limit of the UEG. However, it is not trivial to explain the bad performance of SCAN compared to PBE in the case of a weak perturbation and $q>2q_F$ due to the large number of exact conditions enforced. A more in-depth analysis is needed to explain this behavior.

The presented data along with our assessment constitute an indispensable guide on the choice of the XC functional for KS-DFT calculations when an inhomogeneous electronic structure of WDM is investigated. This is of paramount importance for the diagnostics of XRTS experiments. We highlight the importance of this application by pointing out that KS-DFT results are used to extract electronic parameters like temperature and density from experimental observations.
Besides that, our findings advance our understanding on how well KS-DFT is capable of capturing the electronic structure under WDM conditions. They also point to the parameter space where XC functionals ought to be improved for their use in the WDM application domain.
We highlight the need for XC functionals that are accurate when perturbed electronic states are present. This goes beyond the inclusion of explicit temperature effects in the XC free energy.
This particular outcome of our assessment is valuable for DFT developers.

In the light of the vast amount of available XC functionals, a comprehensive assessment was beyond the scope of this work.
An extensive comparison of available XC functionals shall be presented elsewhere. Instead, in this work, we focus on a representative set of XC functionals $-$ the basic LDA and its common generalizations. Furthermore, we have set up a numerical workflow for benchmarking XC functionals under perturbed electronic structures in WDM.
The presented KS-DFT data, input scripts, and QMC data will be made accessible online~\cite{data}. In doing so, we provide tools to test any existing or newly developed XC functional under perturbed electronic states in WDM. We believe that our work is a valuable contribution that facilitates both the rigorous assessment and construction of XC functionals adapted to the needs of WDM.

\section*{Acknowledgments}
We acknowledge helpful feedback from M. Bussmann. ZM gratefully acknowledges stimulating discussions with Timothy Callow.
This work was funded by the Center for Advanced Systems Understanding (CASUS) which is financed by the German Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science, Art, and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament. We gratefully acknowledge computation time at the Norddeutscher Verbund f\"ur Hoch- und H\"ochstleistungsrechnen (HLRN) under grant shp00026, and on the Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universit\"at Dresden.

\section*{Appendix: KS-DFT Simulation Parameters}\label{sec:app}
All KS-DFT calculations were performed with the GPAW code~\cite{GPAW1, GPAW2, ase-paper, ase-paper2}.
A \emph{k}-point grid of $N_k\times N_k\times N_k$ with $N_k=12$ at $r_s=2$ and $N_k=8$ at $r_s=6$ with a Monkhorst-Pack sampling of the Brillouin zone ($\vec k$-points) was used.
At $T=T_F$, 180 orbitals (with the smallest occupation number of about $10^{-4}$) were used for a total of $N=14$ electrons. The grid spacing was set to $0.15~{\rm \AA}$ for $0.02\leq A\leq 1$ and $r_s=2$, $0.05~{\rm \AA}$  for $A=5$ and $r_s=2$, $0.25~{\rm \AA}$ for $r_s=6$ and $0.05\leq A\leq 0.1$. At $A=5$, $240$ bands were used. Convergence criteria used for the self-consistency cycle: the energy change (last 3 iterations) must be less than $0.5~ {\rm meV}$ per valence electron, the change in integrated absolute value of density change must be less than 0.0001 electrons per valence electron, and the integrated value of the square of the residuals of the Kohn-Sham equations should be less than $4\times 10^{-8}~{\rm eV}^2$ per valence electron (see \href{https://wiki.fysik.dtu.dk/gpaw/documentation/manual.html#manual-convergence}{GPAW documentation}).

In Fig. \ref{fig:Nband_test} we compare results computed using $180$ bands and $280$ bands using LDA and SCAN. From Fig. \ref{fig:Nband_test} we see that there is no notable difference in the results for densities.

In Fig. \ref{fig:NullXC}, we present results for the total density computed using null XC functional in comparison with the data obtained using QMC simulation and  KS-DFT simulation with LDA XC functional. From Fig. \ref{fig:NullXC} we see that it is essential to use XC functional to correctly compute the total density using KS-DFT method.

\begin{figure}
\center
\includegraphics{Nband_test.jpg}
\caption{\label{fig:Nband_test}
Relative deviation in the density $\Delta n/{\rm max }~\delta n  \times 100~\%$ between the KS-DFT data and the reference QMC data for $r_s=2$ with $A=0.1$. The wave numbers of the harmonic perturbation in Eq.~(\ref{eq:H}) are $q_{\rm min}$ (left) and $3 q_{\rm min}$ (right). Top panel: the data computed using LDA XC functional with $180$ and $280$ bands. Bottom panel: the data calculated using SCAN XC functional with $180$ and $280$ bands.}
\end{figure}

\begin{figure}
\center
\includegraphics{Null_xc_A1.jpg}
\caption{\label{fig:NullXC}
The electronic density distribution along the perturbation direction for two different amplitudes $A$ at wave number $q_{\rm min}$ for $r_s=2$ and $T=T_F$.
The comparison of the data computed using LDA XC functional and null XC functional with QMC data shows the importance of the XC functional for the calculation of the electronic density.
}
\end{figure}

\section*{Appendix: PIMC simulation details}\label{sec:app2}

The basic equation behind the PIMC method~\cite{cep,Berne_JCP_1982,Takahashi_Imada_PIMC_1984} is the canonical partition function
\begin{eqnarray}\label{eq:Z}
Z_{\beta,N,V} &=& \frac{1}{N^\uparrow! N^\downarrow!} \sum_{\sigma^\uparrow\in S_N} \sum_{\sigma^\downarrow\in S_N} \textnormal{sgn}(\sigma^\uparrow,\sigma^\downarrow)\\\nonumber & & \times \int d\mathbf{R} \bra{\mathbf{R}} e^{-\beta\hat H} \ket{\hat{\pi}_{\sigma^\uparrow}\hat{\pi}_{\sigma^\downarrow}\mathbf{R}}\ ,
\end{eqnarray}
where $\mathbf{R}=(\mathbf{r}_1,\dots,\mathbf{r}_N)^T$ contains the coordinates of all $N$ particles,  and $\hat{\pi}_{\sigma^\uparrow}$ ($\hat{\pi}_{\sigma^\downarrow}$) is the permutation operator corresponding to a particular element $\sigma^\uparrow$ ($\sigma^\downarrow$) from the permutation group $S_N$. In addition, $\textnormal{sgn}(\sigma^\uparrow,\sigma^\downarrow)$ is the sign function, which is positive (negative) for an even (odd) number of pair exchanges~\cite{Dornheim_permutation_cycles}.
Unfortunately, the matrix elements of the density operator $\hat\rho=e^{-\beta\hat{H}}$ cannot be evaluated in a straightforward way as the kinetic ($\hat{K}$) and potential ($\hat{V}$) contributions to the full Hamiltonian do not commute,
\begin{eqnarray}
\bra{\mathbf{R}}e^{-\beta\hat{H}}\ket{\mathbf{R}} \neq \bra{\mathbf{R}} e^{-\beta\hat{K}} e^{-\beta\hat{V}} \ket{\mathbf{R}}\ .
\end{eqnarray}
As a practical workaround, one uses the exact semi-group property of $\hat{\rho}$, which allows one to re-write Eq.~(\ref{eq:Z}) as
\begin{eqnarray}\label{eq:Z_modified}
Z_{\beta,N,V} &=& \frac{1}{N^\uparrow! N^\downarrow!} \sum_{\sigma^\uparrow\in S_N} \sum_{\sigma^\downarrow\in S_N} \textnormal{sgn}(\sigma^\uparrow,\sigma^\downarrow)\\\nonumber & & \times \int d\mathbf{R}_0\dots d\mathbf{R}_{P-1}
\bra{\mathbf{R}_0}e^{-\epsilon\hat H}\ket{\mathbf{R}_0}\\\nonumber & & \times \bra{\mathbf{R}_1}e^{-\epsilon\hat H}\ket{\mathbf{R}_1} \dots
\bra{\mathbf{R}_{P-1}} e^{-\epsilon\hat H} \ket{\hat{\pi}_{\sigma^\uparrow}\hat{\pi}_{\sigma^\downarrow}\mathbf{R}_0}\ ,
\end{eqnarray}
where we have defined $\epsilon=\beta/P$.
The comparison between Eqs.~(\ref{eq:Z}) and (\ref{eq:Z_modified}) reveals that the partition function has been transformed into a high-dimensional integration over $P$ density matrices that have to be evaluated at $P$ times the original temperature $T$. For a sufficiently large $P$, one can introduce a suitable high-temperature approximation like the simple \emph{primitive factorization} $e^{-\epsilon\hat{H}}\approx e^{-\epsilon\hat{V}}e^{-\epsilon\hat{K}}$. In fact, the associated factorization error decays as $P^{-2}$~\cite{brualla_JCP_2004,sakkos_JCP_2009}, and the convergence in the limit of $P\to\infty$ is ensured by the well-known Trotter formula~\cite{deRaedt_Trotter}. For completeness, we mention that higher-order factorizations of $\hat{\rho}$ are frequently employed in PIMC simulations~\cite{sakkos_JCP_2009,dornheim_njp15,dornheim_cpp_19}, although we do not find them necessary for the present study. In practice, we use $P=200$ primitive high-temperature factors such that any factorization error is substantially below the given level of statistical uncertainty; see the Supplementary Material of Ref.~\onlinecite{Dornheim_PRL_2020} for corresponding convergence plots.

The high dimensionality of Eq.~(\ref{eq:Z_modified}) requires a stochastic evaluation, which can be done efficiently using methods that are based on the celebrated Metropolis algorithm~\cite{metropolis}. Specifically, we employ a canonical adaption~\cite{mezza} of the worm algorithm by Boninsegni and co-workers~\cite{boninsegni1,boninsegni2}. Additional care has to be taken due to the antisymmetric nature of the fermionic density matrix, which can result in positive and negative contributions to Eq.~(\ref{eq:Z_modified}). The corresponding cancellation of positive and negative terms is the origin of the notorious fermion sign problem~\cite{dornheim_sign_problem,troyer}, which leads to an exponential increase in compute time with increasing the system size $N$ or decreasing the temperature $T$.

In practice, this bottleneck is often avoided by imposing the \emph{fixed-node approximation}~\cite{Ceperley1991}. While this allows one to formally remove the sign problem, this advantage comes at the cost of an uncontrolled approximation~\cite{schoof_prl15}. Obviously, this is a disadvantage that severely limits the value of such data to benchmark other approximations like thermal DFT. Therefore, we do not impose any nodal restrictions in our simulations, and our PIMC results are exact within the given level of uncertainty. This is possible for the present parameters by using modern supercomputer clusters, and we have spend $\mathcal{O}(10^5)$ CPUh for the most costly data points with $N=34$ at $r_s=2$ and $\theta=1$.

\section*{Data Availability}
The data supporting the findings of this study are available on the Rossendorf Data Repository (RODARE)~\cite{data}.


\section{Introduction}

\IEEEPARstart{T}{he} Omnidirectional Spherical Robots are the class of mobile robots which are generally having Spherical shell and internal driving components that provide torques required for their rolling motion. Due to a perfectly symmetric geometric structure, SRs acquire multiple advantages over other mobile robots, such as improved mobility, maneuverability, omnidirectionality and the compact structure. The literature survey \cite{crossley2006literature} illustrated the different designs of SRs. Since the first contact wheeled driven SR was developed by Halme et al. in 1996 \cite{halme1996spherical}, the improved two contact wheeled configuration was given by K. Husoy \cite{Husoy2003instrumentation}. A. Bicchi et al. proposed a small resting car-propelled design of SR \cite{bicchi1997introducing}. S. Bhattacharya developed a design that involved a set of two mutually perpendicular rotors attached to the inside of the sphere \cite{bhattacharya2000design}. Joshi and Banavar used four rotors to achieve omnidirectionality in SRs \cite{joshi2010design}. While, Gajamohan et al. \cite{gajamohan2013cubli} used three rotors in a cube robot to rotate and balance omnidirectionally. SR called Rotundus produces pendulum based propulsion \cite{Rotundus}. The Spherobot developed by R. Mukherjee et al. used fixed radially distributed weights along spokes for the actuation \cite{mukherjee1999simple}. The four axle mounted tetrahedral stepper motor driven design was proposed by A. Javadi and P. Mojabi \cite{mojabi2002introducing}, Sang et al. \cite{sang2010modeling} and Tomik et al. \cite{tomik2012design}. Sugiyama et al. proposed a shape-memory alloy to deform and roll wheels and spheres \cite{sugiyama2005circular}, while Wait et al. \cite{wait2010self} and Artusi et al. \cite{artusi2011electroactive}  deform panels on a robot’s shell via air bladders and dielectric actuators.

\begin{figure}[t]
  \begin{center}
  \includegraphics[width=0.70\linewidth]{Images/SPIROS.JPEG}
  \caption{CAD model of the SPIROS, the proposed Omnidirectional Spherical Robot.}\label{SpirosCAD}
  \end{center}
\end{figure}

Recent researches have solved the complexity of the dynamics and control problems of the SRs and found the suitable applications of SRs in surveillance \cite{surveillance}, environmental monitoring \cite{monitoring}, exploration operations in unknown environments \cite{uenvironment}, agriculture \cite{agriculture}. Nevertheless, the vertical motion and the steep inclination maneuver has been an unsolved problem with the existing contact wheeled based designs. Also, the contact wheel's skidding, tipping over, falling or friction with the surface makes the SRs vulnerable or inefficient \cite{friction}. The deformable robot designs does not meet the level of traditional spherical rolling robots, they are somewhat limited by their tethered operation and high voltage requirements. As a result, the fully efficient functionality of traditional SR characteristics cannot be achieved. The design concept to have the sphere split apart to deploy a leg mechanism can't actuate the robot on the vertical surfaces.

Due to NASA’s efforts to create low-cost Mars exploration rovers, there has been some work on wind-powered SRs based on the biomimetic design which is inspired by tumbleweeds \cite{antol2003low} shown in Fig. \ref{wind_spheres}. Such a design would make the robot to go over the larger obstacles than a car-like rover of its size. Other designs have mixed the pendulum design with wind-powered robots, allowing them to steer and stop by using a pendulum that can be retracted into the center of the sphere to “coast” \cite{jones1998beach}. Since then the possibilities of wind-powered SRs have not been fully explored. Air-based propulsion would allow the robot to maneuver over steep as well as vertical surfaces. Moreover, by introducing the air propulsion based drive in SRs, there will be no frictional losses between the driving mechanism and the outer shell, while maintaining a smaller size than the wind spheres, the improved functionality can be achieved.

\begin{figure}[h]
  \begin{center}
  \includegraphics[width=0.9\linewidth]{Images/wind_spheres.PNG}
  \caption{Tumbleweeds inspired wind-powered Wind Spheres}\label{wind_spheres}
  \end{center}
\end{figure}

This paper presents a novel design of SR called SPIROS: The Six Propeller and Intermeshed Rotor based Omnidirectional Spherical Robot. SPIROS uses six intermeshed rotary air thrusters \cite{wildhaber1974rotor} placed in the octahedral arrangement for propulsion. It is omnidirectional (which is able to roll in any direction instantaneously) unlike most of the existing SR designs. The Goldberg Polyhedral \cite{hart2013goldberg} shaped gridshell is used to achieve the outer geometric spherical shape and to ensure the free airflow needed for the thrusters. We have described the mechanical design of the intermeshed rotary thrusters, the outer gridshell and the full robot assembly in detail. Also, as per the proposed classification and the constraints of SRs in \cite{moazami2019kinematics}, the dynamics equations for SRs rolling over 3D terrains are derived. A modified pure pursuit method \cite{corke2017robotics} is utilized for the path tracking control problem of derived models.

The remainder of this paper is organized as follows: Section II elaborates on the working principle of the robot. Section III is dedicated towards the mechanical and electrical design specifications of the robot assembly. In Section IV, the dynamic equations of SPIROS are derived based on the Continuous Rolling Spherical Robots (CR-SR) model, with triple axes rolling (3R-SR), dual axes rolling (2R-SR), and rolling and turning (RT-SR) spherical robots as subcategories. In Section V, a 3D path tracking control scheme has been presented. The simulation results in MATLAB and Simulink are presented in Section VI, followed by the conclusion and future work.

\section{Working Principle}

The propulsion mechanism of SPIROS is based on the principle of the Force Couple or the Pure Moment. According to the Varignon's Second Moment Theorem \cite{VaringnonTheorem}, the torque (moment)  of a couple is independent of the reference point. As shown in Fig.~\ref{tnetlnet}(a), a couple of two skewed clockwise (CW) and counterclockwise (CCW) rotating propellers are placed with $120^{0}$ angles between their axes. The direction of thrust produced by each propeller with equal magnitude $F$ is kept in such a way that the resultant angular momentum $\Vec{L}_{net}$ is neutralized (Fig.~\ref{tnetlnet}(b)), and Force Couple is formed with net torque (moment) $\Vec{\tau}_{net}= \sqrt{3}.l.F$ as shown in Fig.~\ref{tnetlnet}(c).

\begin{figure}[h]
  \begin{center}
  \includegraphics[width=\linewidth]{Diagrams/Torque_angular_momentum.pdf}
  \caption{(a) Arrangement of skewed CW and CCW rotating propeller pairs in Spherical Shell of radius $R$, with the directions of produced thrust $\Vec{F}$, Angular momentum $\Vec{L}$, and the angular rotation $\Vec{\omega}$. (b) Neutralised net angular momentum $\Vec{L}_{net} = 0$. (c) Force Couple with net force $\Vec{F}_{net} = 0$ and torque $\Vec{\tau}_{net} = \sqrt{3}.l.F$.}\label{tnetlnet}
  \end{center}
\end{figure}

As shown in Fig.~\ref{wprinciple}(a), for the controlled maneuver over the slant surface, static friction $\Vec{fr}$ is utilized as an external force acting on the point of contact $P_{0}$. Given a surface slope $\theta$ relative to the gravity vector $\Vec{g}$, mass of the Robot $M$, the torque produced by the thrusters $\Vec{\tau}_{net}$, the angular acceleration of $\ddot{\zeta}$ can be expressed as:

\begin{equation}\label{zetaddot}
    \ddot{\zeta} = -\frac{\tau_{net}+ R.M.g.\sin(\theta)}{I_{0}+M.R^{2}}
\end{equation}

Whereas, $R$ is the radius of the sphere, and $I_{0}$ is the moment of inertia calculated about the axis of the rotation of the robot.

\begin{figure}[h]
  \begin{center}
  \includegraphics[width=\linewidth]{Diagrams/Working_principle.pdf}
  \caption{Schematic diagrams of (a) The basic forces and moments acting on the rolling Robot. (b) Reference frames and generalized robot dynamics parameters.}\label{wprinciple}
  \end{center}
\end{figure}

 Let, at an instant, the robot is rolling over the slant surface heading in $\hat{p}$ direction, $P_{0}$ is an instantaneous point of contact with the surface, $\hat{e}$ be the vector in the surface plane perpendicular to the direction of motion and $\hat{n}$ be the normal vector to the surface from $P_{0}$. Let the moments $\tau_{\hat{p}},\tau_{\hat{e}},\tau_{\hat{n}}$ are applied about axes $\hat{p}$, $\hat{e}$ and $\hat{n}$ to the robot the resultant angular acceleration vector $[\ddot{\zeta}_{\hat{p}},\ddot{\zeta}_{\hat{e}},\ddot{\zeta}_{\hat{n}}]^T $ can be obtained by substituting corresponding 3D quantities in eq. (\ref{zetaddot}), given by;

 \begin{equation}\label{zetaddot3D}
\begin{bmatrix}
\ddot{\zeta}_{\hat{p}}\\
\ddot{\zeta}_{\hat{e}}\\
\ddot{\zeta}_{\hat{n}}
\end{bmatrix}=-(\Tilde{I}_{0}+\Tilde{I}_{R})^{-1}\begin{bmatrix}
\tau_{\hat{p}}\\
\tau_{\hat{e}}-W_{||}R\\
\tau_{\hat{n}}
\end{bmatrix}
\end{equation}

Where, $W_{||}$ is the parallel component of the Weight vector $\Vec{W}$ along $\hat{p}$. The eq. (\ref{zetaddot3D}) is transformed in the terms of three moments $\tau_{\hat{X}}$,$\tau_{\hat{Y}}$ and $\tau_{\hat{Z}}$ along the three orthogonal axes of the robot frame $\mathscr{F}=\{\hat{X},\hat{Y},\hat{Z}\}$ (Fig.~\ref{tnetlnet}(c)). The instantaneous pose of the robot be $\alpha$ (yaw), $\beta$ (pitch) and $\gamma$ (roll). As shown in Fig.~\ref{wprinciple}(b), The resulting angular accelerations about arbitrary three mutually perpendicular axes of the frame $\mathscr{T}=\{\hat{h},\hat{l},\hat{n}\}$ are $\ddot{\Phi}$ along longitudinal axis $\hat{h}$, $\ddot{\Theta}$ along lateral axis $\hat{l}$, and $\ddot{\Psi}$ along normal axis $\hat{n}$ of the robot. The eq. (\ref{zetaddot3D}) is transformed to obtain angular acceleration vector $[\ddot{\Phi},\ddot{\Theta},\ddot{\Psi}]^{T}$ is given by:

\begin{equation}\label{genangddot}
    \begin{aligned}
    \begin{bmatrix}
     \ddot{\Theta}\\
     \ddot{\Phi}\\
     \ddot{\Psi}
    \end{bmatrix} = -R_{\xi\hat{n}}^{T}\Tilde{I}^{-1}R_{\hat{p}}R_{\alpha\beta\gamma}^{T}\begin{bmatrix}
     \tau_{\hat{X}}\\
     \tau_{\hat{Y}}\\
     \tau_{\hat{Z}}
    \end{bmatrix}-R_{\xi\hat{n}}^{T}\Tilde{I}^{-1}\begin{bmatrix}
     0\\
     W_{||}R\\
     0
    \end{bmatrix}
    \end{aligned}
\end{equation}

Where $R_{\alpha\beta\gamma}$ (eq. (\ref{Rabg})) is a general rotation matrix (Appendix \ref{'Rabc'}) gives rotational transform from robot frame $\mathscr{F}$ to the world frame $\mathscr{W}$. The matrices $R_{\xi\hat{n}}$ (eq. (\ref{Rxin})) and $R_{\hat{p}}$ (eq. (\ref{pen}) ) are the rotational transforms from the frame $\mathscr{T}$ to the frame $\{\hat{p},\hat{e},\hat{n}\}$ and the frame $\{\hat{p},\hat{e},\hat{n}\}$ to the world frame $\mathscr{W}$ respectively. (Appendix \ref{'ddot'}). $\Tilde{I} = \Tilde{I}_{0}+\Tilde{I}_{R})$ is the Inertia Tensor matrices calculated about $P_{0}$ illustrated in Appendix \ref{'Inot'}.

\section{Hardware Analysis}

To conceptualize the basic dynamic equation (\ref{genangddot}) of the robot with the required configurations of moments and thrusts discussed in Section II, the mechanical arrangement of the SPIROS is designed as follows.

\subsection{Intermeshing Rotary Thruster}

The propeller thrust is the main driving force of the robot. The clockwise and anticlockwise rotation of the propeller pair shown in Fig~\ref{tnetlnet}, is incorporated by using the bevel gear \cite{reed1908bevel} arrangement. Three bevel gears with equal gear ratios are required with an angle $60^0$ with each other to transfer the driving rotational motion into the two $120^0$ skewed anti-rotating shafts. The schematic of the gearbox and the axle support with 3mm bearings are shown in the Fig.~\ref{Thruster}(b). The required outer shell diameter is minimized by the skewed axis of the shafts.

The thruster is actuated with the high-speed motors such as BLDC motors with its suitable ESC (Electronic Speed Controller) \cite{xia2012permanent}. The mechanical mounts for the driving motor and its speed controller are as shown in the Fig.~\ref{Thruster}(a). To obtain the required quantity of air thrust, the 10" four blades propellers are used (Shown in Fig.~\ref{Thruster}(c)). The overall structure of the Intermeshing Rotary Thruster is   modularly designed (as shown in Fig.~\ref{Thruster}) to repetitively construct and arrange in six-thruster configuration shown in Fig.~\ref{OverallConstruction}(a).

\begin{figure}[h]
  \begin{center}
  \includegraphics[width=0.95\linewidth]{Diagrams/Thruster_assembly.pdf}
  \caption{Schematic diagram of Intermeshing Rotary Thruster, (a) The BLDC motor with its Electronic Speed Controller (ESC). (b) The Bevel Gear mechanism. (c) 10" 4-blade Propeller}\label{Thruster}
  \end{center}
\end{figure}

\subsection{Inner Support Cube}

The inner support cube shown in Fig.~\ref{InnerCube} functions as the main supporting frame of the Robot's internal components. It comprises of thruster mounts, outer shell supports and a central hollow frame. The central hollow space provides room for the mechanical components as well as electronics other than thrusters.

Due to the central position of the structure, the Robot's non-inertial reference frame  $\mathscr{F}$ is defined along the three orthogonal vectors $\hat{X}$,$\hat{Y}$ and $\hat{Z}$ passing through center of the cube $O$ and centers of the respective faces.

\begin{figure}[h]
  \begin{center}
  \includegraphics[width=0.5\linewidth]{Diagrams/Inner_Cube.pdf}
  \caption{Schematic diagram of Inner Support Cube comprising of (a) The Thruster Mounts, (b) The Internal Hollow Frame, and (c) The Outer shell supports.}\label{InnerCube}
  \end{center}
\end{figure}

The outer gridshell (Section III C.) is rigidly attached to the shell supports of the central cube. Due to the overall rolling motion of the robot, the internal framework, as well as whole thruster assembly, rolls with it.  Therefore, the transformation between the world frame $\mathscr{W}=\{\hat{x},\hat{y},\hat{z}\}$ and robot frame $\mathscr{F}$ is given by:

\begin{figure*}[!b]
  \begin{center}
  \includegraphics[width=\linewidth]{Diagrams/Overall_Robot_Assembly.pdf}
  \caption{(a) Diametric View of the internal arrangement of the SPIROS, Schematics of the overall assemblies of the SPIROS in (b) Isometric View and (c) Front View.}\label{OverallConstruction}
  \end{center}
\end{figure*}

\begin{equation}
     \begin{bmatrix}
     \mathscr{W}_{3\times1}\\
     1
     \end{bmatrix}=\begin{bmatrix}
     R_{\alpha\beta\gamma} & P_{0}\\
     0&1
     \end{bmatrix}.\begin{bmatrix}
     \mathscr{F}_{3\times1}\\
     1
     \end{bmatrix}
 \end{equation}

\subsection{Outer Spherical Gridshell}

To ensure the free air circulation within the thrusters as well as to maintain structural rigidity and defined spherical geometry, the sphere-shaped perforated grid is used as an outer shell of the robot. Taking into account the structural properties of hexagon-pentagonal spherical shapes discussed in \cite{malek2013structural}\cite{martin1975ellipse}, such a grid shape offers rigidity and the free passage of the gases.

\begin{figure}[h]
  \begin{center}
  \includegraphics[width=\linewidth]{Diagrams/Outer_Shell.pdf}
  \caption{(a) The outer gridshell of the SPIROS, (b) Steps of building Goldberg Polyhedral grid.}\label{OuterShell}
  \end{center}
\end{figure}

The required shape of the grid considerably resembles the Goldberg polyhedron \cite{hart2013goldberg}. The method suggested in \cite{1937104} is applied to design the outer gridshell as illustrated in Fig.~\ref{OuterShell}. To get sufficient friction for pure rolling, the outer lining of the shell is coated with high friction coefficient materials as given in \cite{majidi2006high}.

\subsection{Electrical Architecture}

The electrical architecture illustrated in Fig.~\ref{ElectArch} is utilized in the electrical system of the SPIROS. The reference path acquired from the user input is interpreted and the required torque vector is calculated. The angular speed signals are given to the thruster pairs via motor drivers. The IMU sensor \cite{chen1994gyroscope} and the Hall effect sensors are used to receive the robot's instantaneous pose $\{\alpha,\beta,\gamma\}$ and the output thrust magnitude's feedback during the course of motion.

\begin{figure}[ht]
  \begin{center}
  \includegraphics[width=\linewidth]{Diagrams/Electrical_Architecture.pdf}
  \caption{Electrical Architecture of the SPIROS.}\label{ElectArch}
  \end{center}
\end{figure}

\subsection{Overall Mechanical Assembly}

The thrusters (Section III A.) are attached to the central cubical frame (Section III B.) in the configuration schematically shown in Fig.~\ref{OverallConstruction}(a). To ensure the dynamic conditions proposed in the working principle (Section II.) of the robot, the magnitude of the thrust generated by the diametrically opposite thrusters are kept equal by the speed of the respective actuators. This arrangement enables the generation of the controlled moments $\tau_{\hat{X}}$, $\tau_{\hat{Y}}$ and $\tau_{\hat{Z}}$ along the axes of the robot frame $\mathscr{F}=\{\hat{X},\hat{Y},\hat{Z}\}$.  Table.~\ref{table:ThrustMoments} gives the geometric coordinates, thrust generated and the resultant torques by each thruster with respect to the robot frame $\mathscr{F}$.

\begin{table}[h!]
\centering
\caption{Geometric positions, the thrust generated and the resultant torques by each thruster with respect to the robot frame $\mathscr{F}$.}
\begin{tabular}{ c c c c }
\hline
No. & Position $\Vec{R}$ & Thrust $\Vec{F}$& Torque $\tau=\Vec{R}\times\Vec{F}$\\
\hline
\hline

1. & $[0,l,0]^{T}$ & $[0,0,F_{\hat{X}}]^{T}$ & $[lF_{\hat{X}},0,0]^{T}$\\
2. & $[0,-l,0]^{T}$ & $[0,0,-F_{\hat{X}}]^{T}$ & $[lF_{\hat{X}},0,0]^{T}$\\
3. & $[0,0,l]^{T}$ & $[F_{\hat{Y}},0,0]^{T}$ & $[0,lF_{\hat{Y}},0]^{T}$\\
4. & $[0,0,-l]^{T}$ & $[-F_{\hat{Y}},0,0]^{T}$ & $[0,lF_{\hat{Y}},0]^{T}$\\
5. & $[l,0,0]^{T}$ & $[0,F_{\hat{Z}},0]^{T}$ & $[0,0,lF_{\hat{Y}}]^{T}$\\
6. & $[-l,0,0]^{T}$ & $[0,-F_{\hat{Z}},0]^{T}$ & $[0,0,lF_{\hat{Y}}]^{T}$\\
\hline
\end{tabular}
\label{table:ThrustMoments}
\end{table}

\section{Kinematics and Dynamics Modelling}

Consider a spherical SPIROS robot rolling over a 3D surface $\sigma$ (Appendix \ref{'ddot'}), as shown in Fig.~\ref{Dynamics}. To derive the kinematics and the dynamics equations of rolling motion, the reference frames and the points are defined as follows. The position of the contact point of the robot with $\sigma$ is represented as $P_{0}$, which instantaneously indicates the position coordinates $[x(t),y(t),z(t)]^T$ of the robot in the inertial world frame $\mathscr{W}=\{\hat{x},\hat{y},\hat{z}\}$. The moments generated by the thruster assembly ($[\tau_{\hat{X}},\tau_{\hat{Y}},\tau_{\hat{Z}}]^T$) are considered along the axes of the non-inertial robot frame $\mathscr{F}=\{\hat{X},\hat{Y},\hat{Z}\}$. The robot frame rotates about the sphere center $O$ as the robot performs the rolling motion.

\begin{figure}[h!]
  \begin{center}
  \includegraphics[width=0.6\linewidth]{Diagrams/dynamics.pdf}
  \caption{Schematic diagram of reference frames and points that are used to derive the kinematics and dynamics equation of the rolling SPIROS on a 3D surface $\sigma$.}\label{Dynamics}
  \end{center}
\end{figure}

The frame $\mathscr{T}=\{\hat{h},\hat{l},\hat{n}\}$ is defined with its axes $\hat{h}$ and $\hat{l}$ laying on the imaginary tangent plane $S_{\sigma}$ passing through $P_{0}$. The normal axis $\hat{n}$ of the robot is defined perpendicular to the  $S_{\sigma}$ coinciding with the center of the sphere $O$, while lateral and longitudinal axes are defined parallel to the $S_{\sigma}$, That is, they do not rotate through rotation of the sphere about axes parallel to the $S_{\sigma}$, but they rotate accordingly when the sphere rotates about the normal axis $\hat{n}$. The rolling angle $\Theta$, tilting angle $\Phi$ and the turning angle $\Psi$ are as per defined in Section II.

Assumptions: The mathematical model of the SPIROS is derived considering the following assumptions:

\begin{enumerate}
    \item The surface $\sigma$ is rough enough to allow the sphere to roll over without slipping, skidding or falling; That is, without losing the traction with the surface.
    \item The surface $\sigma:z=f(x,y)$ is a continuous function of class $C_{2}$ or higher \cite{barsky1989geometric}.
    \item The outer gridshell is rigid, and the radius of the curvature of the surface $\sigma$ is never smaller than the radius of the gridshell, therefore, the robot remains in the single point of contact with the surface.
    \item The moments generated and the frictional force are always enough to provide the required driving force to the robot.
\end{enumerate}

Regarding to the kinematics equations derived by Moazami et al. \cite{moazami2019kinematics} based on the different categories of the SRs rolling over the 3D terrains, the transformation between the world frame $\mathscr{W}$ co-ordinates and the tangent plane frame $\mathscr{T}$ co-ordinates is given by:

\begin{equation}
    \begin{bmatrix}
    \mathscr{W}_{3\times1}\\
    1
    \end{bmatrix}=\begin{bmatrix}
    a_{11}&a_{12}&a_{13}&x\\
    a_{21}&a_{22}&a_{23}&y\\
    a_{31}&a_{32}&a_{33}&z\\
    0&0&0&1\\
    \end{bmatrix}\begin{bmatrix}
    \mathscr{T}_{3\times1}\\
    1
    \end{bmatrix}
\end{equation}

and,

\begin{equation}\label{hlnTrans}
    [\hat{h},\hat{l},\hat{n}]=\begin{bmatrix}
    a_{11}&a_{12}&a_{13}\\
    a_{21}&a_{22}&a_{23}\\
    a_{31}&a_{32}&a_{33}\\
    \end{bmatrix}
\end{equation}

where the elements $a_{jk}$ are given by:

\begin{equation}
\begin{aligned}
\begin{array}{l}
a_{11}=(1-s^{2}f_{x}^{2}(1-s_{n}))C\Psi+s^{2}f_{x}f_{y}(1-s_{n})S\Psi\\
a_{12}=(1-s^{2}f_{x}^{2}(1-s_{n}))S\Psi-s^{2}f_{x}f_{y}(1-s_{n})C\Psi\\
a_{13}=-s_{n}f_{x}\\
a_{21}=-s^{2}f_{x}f_{y}(1-s_{n})C\Psi-(1-s^{2}f_{y}^{2}(1-s_{n}))S\Psi \\
a_{22}=-s^{2}f_{x}f_{y}(1-s_{n})S\Psi+(1-s^{2} f_{y}^{2}(1-s_{n}))C\Psi\\
a_{23}=-s_{n}f_{y}\\
a_{31}=s_{n}f_{x}C\Psi-s_{n}f_{y}S\Psi\\
a_{32}=s_{n}f_{x}S\Psi+s_{n}f_{y}C\Psi\\
a_{33}=s_{n}\\
\end{array}
\end{aligned}
\end{equation}

while, $C\Psi$,~$S\Psi$,~$f_x$,~$f_y$,~$s$ and $s_{n}$ are the terms defined in Appendix \ref{'ddot'}.

The dynamic model of the SPIROS is derived considering the kinematic models of the three sub-categories of the Continuous Rolling Spherical Robot (CR-SR) Model borrowed from \cite{moazami2019kinematics}.

\subsection{Triple-Axis Rolling Spherical Robot (3R-SR) Model}

This is the most general category of the CR-SR. In this model, the SR is considered to rotate about all of its three axes, which is a case for the SPIROS. The kinematical model for the 3R-SR given in \cite{moazami2019kinematics} is as per the eq. (\ref{3RSRkinematics}).

\begin{equation}\label{3RSRkinematics}
    \begin{bmatrix}
    \Dot{x}\\
    \Dot{y}\\
    \Dot{z}
    \end{bmatrix}=R\mathscr{L}\begin{bmatrix}
    \Dot{\Theta}\\
    \Dot{\Phi}
    \end{bmatrix},\mathscr{L}=R\begin{bmatrix}
    a_{11}&-a_{12}\\
    a_{21}&-a_{22}\\
    a_{31}&-a_{32}
    \end{bmatrix}
\end{equation}

Therefore,

\begin{equation}\label{linv}
    \begin{bmatrix}
    \Dot{\Theta}\\
    \Dot{\Phi}
    \end{bmatrix}=\mathscr{L}^{\dagger}\begin{bmatrix}
    \Dot{x}\\
    \Dot{y}\\
    \Dot{z}
    \end{bmatrix},\mathscr{L}^{\dagger}=\frac{1}{R}\begin{bmatrix}
    \frac{a_{22}}{s_{n}^{2}}&-\frac{a_{12}}{s_{n}^{2}}&0\\
    0&-\frac{a_{31}}{f_{x}s_{n}}&\frac{a_{21}}{f_{x}s_{n}}
    \end{bmatrix}
\end{equation}

By differentiating eq.(\ref{3RSRkinematics}) on both sides and combining it with eq.(\ref{linv}) leads to:

\begin{equation}\label{3RSRdynamics}
   \begin{bmatrix}
    \Ddot{x}\\
    \Ddot{y}\\
    \Ddot{z}
    \end{bmatrix}=\Dot{\mathscr{L}}\mathscr{L}^{\dagger}\begin{bmatrix}
    \Dot{x}\\
    \Dot{y}\\
    \Dot{z}
    \end{bmatrix}+\begin{bmatrix}
    \mathscr{L}&0_{3\times1}
    \end{bmatrix}\begin{bmatrix}
    \Ddot{\Theta}\\
    \Ddot{\Phi}\\
    \Ddot{\Psi}
    \end{bmatrix}
\end{equation}

Therefore, from eq. (\ref{3RSRdynamics}) and eq. (\ref{genangddot}), the dynamics model is written in the following canonical form:

\begin{equation}\label{DynModel1}
    \Ddot{\mathscr{X}}_{\mathscr{W}}
    =\mathscr{A}_{3\times3}\Dot{\mathscr{X}}_{\mathscr{W}}+\mathscr{B}_{3\times3}\tau_{\mathscr{F}} + \mathscr{D}_{3\times1}
\end{equation}

Whereas,

\begin{equation}
    \begin{aligned}
        \begin{array}{l}
        \vspace{2mm}
            \mathscr{X}_{\mathscr{W}}=\begin{bmatrix}
            x\\
            y\\
            z
            \end{bmatrix},\tau_{\mathscr{F}}=\begin{bmatrix}
            \tau_{\hat{X}}\\
            \tau_{\hat{Y}}\\
            \tau_{\hat{Z}}
            \end{bmatrix}\\\vspace{5mm}
             \mathscr{A}_{3\times3} = -\Dot{\Psi}\begin{bmatrix}
             \frac{a_{12}a_{22}}{s_{n}^{2}}&-\frac{a_{11}a_{31}}{f_xs_{n}}-\frac{a_{12}^{2}}{s_{n}^{2}}&\frac{a_{11}a_{21}}{f_xs_{n}}\\
             \frac{a_{22}^{2}}{s_{n}^{2}}&-\frac{a_{21}a_{31}}{f_xs_{n}}-\frac{a_{12}a_{22}}{s_{n}^{2}}&\frac{a_{21}^{2}}{f_xs_{n}}\\
             \frac{a_{22}a_{32}}{s_{n}^{2}}&-\frac{a_{31}^{2}}{f_xs_{n}}-\frac{a_{21}a_{32}}{s_{n}^{2}}&\frac{a_{21}a_{31}}{f_xs_{n}}\\
             \end{bmatrix},\\\vspace{3mm}
             \mathscr{B}_{3\times3}=-\begin{bmatrix}
    \mathscr{L}&0_{3\times1}
    \end{bmatrix}R_{\xi\hat{n}}^{T}\Tilde{I}^{-1}R_{\hat{p}}R_{\alpha\beta\gamma}^{T},\\
             \mathscr{D}_{3\times1}=-\begin{bmatrix}
    \mathscr{L}&0_{3\times1}
    \end{bmatrix}R_{\xi\hat{n}}^{T}\Tilde{I}^{-1}\begin{bmatrix}
    0\\
    \frac{s_nMgR}{s}\\
    0
    \end{bmatrix}
        \end{array}
    \end{aligned}
\end{equation}

\subsection{Dual-Axis Rolling Spherical Robot (2R-SR) Model}

In this model, the robot is considered to be able to roll about its two main axis. As the rotation about the vertical axis is restricted, it is assumed that $\Dot{\Psi}=0$, therefore the eq. (\ref{DynModel1}) can be modified as:

\begin{equation}\label{DynModel2}
    \Ddot{\mathscr{X}}_{\mathscr{W}}
    =\mathscr{B}_{3\times3}\tau_{\mathscr{F}} + \mathscr{D}_{3\times1}
\end{equation}

Hence, the eq. (\ref{DynModel2}) represents the 2R-SR dynamic model of the SPIROS.

\subsection{Rolling and Turning Spherical Robot (RT-SR) Model}

In addition to the rolling action, RT-SRs are supposed to turn about their vertical axis to change their moving direction. The robot mechanism is assumed to provide angular velocities of $\Dot{\Theta}$ and $\Dot{\Psi}$ while $\Dot{\Phi}=0$. Therefore the eq. (\ref{3RSRkinematics}) and (\ref{linv}) can be rewritten as:

\begin{equation}\label{RTSRkinematics}
    \begin{bmatrix}
    \Dot{x}\\
    \Dot{y}\\
    \Dot{z}
    \end{bmatrix}=R\begin{bmatrix}
    a_{11}\\
    a_{21}\\
    a_{31}
    \end{bmatrix}\Dot{\Theta},~\Dot{\Theta}=[a_{11},a_{21},a_{31}]\begin{bmatrix}
    \Dot{x}\\
    \Dot{y}\\
    \Dot{z}
    \end{bmatrix}
\end{equation}

Hence, by differentiating eq. (\ref{RTSRkinematics}) we get;

\begin{equation}\label{RTSTdynamics}
 \begin{bmatrix}
    \Ddot{x}\\
    \Ddot{y}\\
    \Ddot{z}
    \end{bmatrix}=-\Dot{\Psi}\begin{bmatrix}
    a_{12}\\
    a_{22}\\
    a_{32}
    \end{bmatrix}[a_{11},a_{21},a_{31}]\begin{bmatrix}
    \Dot{x}\\
    \Dot{y}\\
    \Dot{z}
    \end{bmatrix}+R\begin{bmatrix}
    a_{11}\\
    a_{21}\\
    a_{31}
    \end{bmatrix}\Ddot{\Theta}
\end{equation}

Therefore, the matrices $\mathscr{A}$, $\mathscr{B}$ and $\mathscr{D}$ of the canonical form (eq. (\ref{DynModel1})) are modified as given below.

\begin{equation}\label{dyn3par}
    \begin{aligned}
        \begin{array}{l}
        \vspace{3mm}
             \mathscr{A}_{3\times3} = -\Dot{\Psi}\begin{bmatrix}
             a_{12}a_{12}&a_{12}a_{21}&a_{12}a_{31}\\
             a_{11}a_{22}&a_{21}a_{22}&a_{22}a_{31}\\
             a_{12}a_{31}&a_{22}a_{31}&a_{31}a_{32}\\
             \end{bmatrix},\\\vspace{3mm}
             \mathscr{B}_{3\times3}=-R.\begin{bmatrix}
    a_{11}&0&0\\
    a_{21}&0&0\\
    a_{31}&0&0
    \end{bmatrix}R_{\xi\hat{n}}^{T}\Tilde{I}^{-1}R_{\hat{p}}R_{\alpha\beta\gamma}^{T},\\
             \mathscr{D}_{3\times1}=-R.\begin{bmatrix}
    a_{11}&0&0\\
    a_{21}&0&0\\
    a_{31}&0&0
    \end{bmatrix}R_{\xi\hat{n}}^{T}\Tilde{I}^{-1}\begin{bmatrix}
    0\\
    \frac{s_nMgR}{s}\\
    0
    \end{bmatrix}
        \end{array}
    \end{aligned}
\end{equation}

\section{Path Tracking Control}

In this section, the pure pursuit method is utilized to design the path tracking controls for the SPIROS based on the three dynamics models derived in Section IV.

Considering $\mathscr{P}=[x_d(t),y_d(t),z_d(t)]^T$ as a desired trajectory laying on $\sigma$, which is assumed to be continuous; that is,$\parallel[\Ddot{x}_d,\Ddot{y}_d,\Ddot{z}_d]^T\parallel$ is bounded for $t>0$. The stable error term $e_t$ with the deviation angle $\zeta$ between the heading direction $\hat{h}$ is defined as:

\begin{equation}
    e_t=[x_d(t),y_d(t),z_d(t)]^T-P_0
\end{equation}

The objective of this algorithm is to converge the
tracking error $e_t$ along with $\zeta$ to zero. The given path-planning expressions computes values for angular accelerations $\Ddot{\Theta}$, $\Ddot{\Phi}$ and $\Ddot{\Psi}$. Further, the eq. (\ref{genangddot}) is used to calculate the respective torques $\tau_{\hat{X}}$,$\tau_{\hat{Y}}$ and $\tau_{\hat{Z}}$.

\subsection{3R-SR Mode Path Tracking Control}

As 3R-SRs dynamic model is capable of providing $\Ddot{\Theta}$ $\Ddot{\Phi}$ and $\Ddot{\Psi}$, the controller is designed to provide the required values accordingly as follows:

\begin{equation}
\begin{array}{l}
\ddot{\Theta}=k_{\Theta} \frac{\|e_{t}\|}{(k_{e}+\|e_{t}\|)} \cos (\zeta)+\frac{1}{R}\|[\ddot{x}_{d}, \ddot{y}_{d}, \ddot{z}_{d}]^{T}\| \\
\ddot{\Phi}=-(k_{\Phi} \frac{\|e_{t}\|}{(k_{e}+\|e_{t}\|)} \sin (\zeta)+k_{\Phi_{2}} \sin (\zeta)) \\
\ddot{\Psi}=k_{\Psi} \zeta
\end{array}
\end{equation}

\subsection{2R-SR Mode Path Tracking Control}

As in 2RSR model, the rotation about the vertical axis is restricted, it can be implied that the robot can only use $\Ddot{\Theta}$ and $\Ddot{\Phi}$ to move towards the target. Therefore, the controller is given by:

\begin{equation}
\begin{array}{l}
\ddot{\Theta}=k_{\Theta_{1}} \frac{\|e_{t}\|}{(k_{e}+\|e_{t}\|)} \cos (\zeta)+k_{\Theta_{2}} \cos (\zeta) \\
\ddot{\Phi}=-(k_{\Phi_{1}} \frac{\|e_{t}\|}{(k_{e}+\|e_{t}\|)} \sin (\zeta)+k_{\Phi_{2}} \sin (\zeta)) \\
\ddot{\Psi}=0
\end{array}
\end{equation}

\begin{figure*}[hb]
  \begin{center}
  \includegraphics[width=0.8\linewidth]{Diagrams/Ramp_simulation.pdf}
  \caption{Robot ascending over the ramp with given initial poses (a)$\{\alpha=0^r,\beta=0^r,\gamma=0^r\}$, (b) $\{\alpha=\frac{\pi}{3}^r,\beta=\frac{\pi}{3}^r,\gamma=\frac{\pi}{3}^r\}$ and (c)$\{\alpha=\frac{\pi}{2}^r,\beta=\frac{\pi}{3}^r,\gamma=\frac{\pi}{6}^r\}$}\label{RampSimulation}
  \includegraphics[width=0.85\linewidth]{Diagrams/Ramp_plots.pdf}
  \caption{Robot's generated moments $\tau_{\hat{X}}$, $\tau_{\hat{Y}}$ and $\tau_{\hat{Z}}$ along the axes $\hat{X}$,$\hat{Y}$ and $\hat{Z}$ of the robot frame $\mathscr{F}$.}\label{RamptorqePlots}
  \end{center}
\end{figure*}

\subsection{RT-SR Mode Path Tracking Control}

As per the RT-SR model the robot can move towards the target point by compensating their deviation error $\zeta$, through turning action and then using $\ddot{\Theta}$ to approach and follow the target. This method is used to design the controller to provide the required angular accelerations as the following:

\begin{equation}
\begin{array}{l}
\ddot{\Theta}=k_{\Theta} \frac{\|e_{t}\|}{(k_{e}+\|e_{t}\|)} \cos(\zeta)+\frac{1}{R}\|[\ddot{x}_{d}, \ddot{y}_{d}, \ddot{z}_{d}]^{T}\|\\
\ddot{\Phi}=0\\
\ddot{\Psi}=k_{\Psi} \zeta
\end{array}
\end{equation}

\section{Simulation Results}

The SPIROS's mathematical model has been simulated in MATLAB and Simulink to verify the dynamics and the controllers presented in Sections IV and V. The physical parameters of the robot are normalized as radius $R = 1m$ and mass $M = 0.5 kg$ (considered as a homogenous sphere) for the simulations. Two operations of the SPIROS are been simulated for the analysis are as follows.

The robot is raised onto the rough ramp (providing enough frictional force for pure rolling) having a slope $\frac{\pi}{8}^{r}$ through the moments generated by thrusters. The experiment is repeated with the different initial poses $\{\alpha,\beta,\gamma\}$ of the robot to analyze the effect on the required torques. The time-lapsed positions of the simulated robot along with its trajectories with initial poses $\{\alpha=0^r,\beta=0^r,\gamma=0^r\}$, $\{\alpha=\frac{\pi}{3}^r,\beta=\frac{\pi}{3}^r,\gamma=\frac{\pi}{3}^r\}$ and $\{\alpha=\frac{\pi}{2}^r,\beta=\frac{\pi}{3}^r,\gamma=\frac{\pi}{6}^r\}$ are shown in Fig.~\ref{RampSimulation} (a),(b) and (c) respectively. Whereas, the corresponding torque magnitudes $\tau_{\hat{X}}$, $\tau_{\hat{Y}}$ and $\tau_{\hat{Z}}$ taken positive along axes $\hat{X}$,$\hat{Y}$, $\hat{Z}$ of robot frame $\mathscr{F}$ are given in the Fig.~\ref{RamptorqePlots} (a),(b) and (c).

Further, the robot is traversed over the continuous rough 3D surface $\sigma$ on the targeted trajectory obtained by projecting a circle over $\sigma$ such that:

\begin{equation}\label{trajectory}
\begin{array}{l}
    x=5+\sin\left(\frac{2\pi}{10}t\right)\\
    y=5+\cos\left(\frac{2\pi}{10}t\right)
\end{array}
\end{equation}

Where the terrain surface, $\sigma$ is defined as the following:

\begin{equation}\label{sigma}
    z=0.5\left(\cos\left(\frac{2\pi}{8}x\right)+\cos\left(\frac{2\pi}{8}y\right)-2\right)
\end{equation}

With controller design parameters illustrated in Section V are tuned by the trial-and-error method, the robot successfully traced the circular trajectory (eq.\ref{trajectory}) shown by 3D time-lapse and trajectory plots in Fig.~\ref{Circle_Simulation}(a). Also, the magnitudes of the torque applied to the robot during the operation are plotted versus time in Fig.~\ref{Circle_Simulation}(b).

\begin{figure*}[ht!]
  \begin{center}
  \includegraphics[width=0.85\linewidth]{Diagrams/Circle_simulation.pdf}
  \caption{(a) Robot traversing over a 3D surface sigma in a smooth circular path. (b) Applied torques $\tau_{\hat{X}}$,$\tau_{\hat{Y}}$,$\tau_{\hat{Z}}$ by the thrusters. }\label{Circle_Simulation}
  \end{center}
\end{figure*}

\section{Conclusion and Futureworks}

In this paper, the proposed SPIROS-Robot mechanism incorporating a novel six thruster propulsion system provides several advantages over the existing SRs. The working principle of the robot is explained. Then the design details of the mechanical structure and the electrical system discussed in detail.  By utilizing the present kinematics models for the SRs, the dynamic equations are derived along with the control schemes based on pure pursuit algorithm are given for the SPIROS. Finally, the derived dynamic equations of the SPIROS are mathematically modeled and numerically simulated in MATLAB and Simulink, showing the expected results.

The future work involves including the effects of the outer shell hindrance and the propeller interaction on the thrusts generated in the mathematical model by the advanced airflow analysis. Also improving on the control scheme to design a robust closed-loop controller. The next step is to develop a hardware prototype of the SPIROS having all the functionalities proposed in this paper.

\appendices

\section{General Rotations}
\label{'Rabc'}

The following matrix $R_{\alpha\beta\gamma}$ represents a rotation whose yaw, pitch, and roll angles are $\alpha$, $\beta$ and $\gamma$ respectively.

\begin{equation}\label{Rabg}
    R_{\alpha\beta\gamma}=\begin{bmatrix}
      C\alpha C\beta & C\alpha C\beta S\gamma-S\alpha C\gamma & C\alpha S\beta C\gamma+S\alpha S\gamma \\
        S\alpha C\beta & S\alpha S\beta S\gamma+C\alpha C\gamma & S\alpha S\beta C\gamma-C\alpha S\gamma \\
      -S\beta & C\beta S\gamma & C\beta C\gamma
    \end{bmatrix}
\end{equation}

Where $S\theta$ and $C\theta$ depicts $\sin(\theta)$ and $\cos(\theta)$ respectively. The pre-multiplication of $R_{\alpha\beta\gamma}$ with the 3D vector will produce the intrinsic rotated vector whose Tait–Bryan angles are $\alpha$, $\beta$ and $\gamma$ about axes $\hat{z}$, $\hat{y}$, $\hat{x}$, respectively.

\section{}
\label{'ddot'}

According to Fig.~\ref{wprinciple}, $\Theta$ is the sphere's rolling angle, that is, the rotation angle of the sphere about its lateral axis. $\Phi$ is the sphere's tilting angle, which is, the rotation angle of the sphere about its longitudinal axis, and $\Psi$ is the sphere's turning angle about the robot-side normal axis. Let $\sigma$ be a $3\mathrm{D}$ surface defined as $\sigma: F(x, y, z)=0$ in the form of $f(x, y)-z=0$ implying that for any $x$ and $y$, the elevation of the surface $\sigma$ along the vertical axis of $\hat{z}$ is $z=f(x, y)$.

The normal vector of surface $\sigma$ at position $P_{0}$ is given by

\begin{equation}\label{n}
\hat{n}=-\frac{\nabla f}{\parallel\nabla f \parallel}=s_{n}[-f_{x},-f_{y},1]^{T}
\end{equation}

where $f_{x}=\frac{\partial f(x, y)}{\partial x}$ and $f_{y}=\frac{\partial f(x, y)}{\partial y}$ are the partial derivative functions and $s_{n}=(f^{2}_{x}+f^{2}_{y}+1)^{-\frac{1}{2}}$. As $\hat{e}$ is perpendicular to the direction of motion and the normal $\hat{n}$.

\begin{equation}\label{e}
  \hat{e}=\hat{n}\times\hat{z} = s[f_{y},-f_{x},0]^{T}
\end{equation}

whereas, $s=(f^{2}_{x}+f^{2}_{y})^{-\frac{1}{2}}$. The vector $\hat{p}$ is along the direction of motion, hence it is perpendicular to the $\hat{e}$ and $\hat{n}$.

\begin{equation}\label{p}
  \hat{p}=\hat{e}\times\hat{n}=ss_{n}[-f_{x},-f_{y},-f^{2}_{x}-f^{2}_{y}]^{T}
\end{equation}

therefore, defining:

\begin{equation}\label{pen}
  R_{\hat{P}}=[\hat{p},\hat{e},\hat{n}]=\begin{bmatrix}
  -ss_{n}f_{x} & sf_{y} & -s_{n}f_{x}\\
  -ss_{n}f_{y} & -sf_{x} & -s_{n}f_{y}\\
  -ss_{n}(f^{2}_{x}+f^{2}_{y}) & 0 & s_{n}
  \end{bmatrix}
\end{equation}

According to the eq. (\ref{hlnTrans}), the vector $\hat{l}$ is given by:

\begin{equation}\label{l}
  \hat{l}=\begin{bmatrix}
  (1-s^{2}f^{2}_{x}(1-s_{n}))S\Psi-s^{2}f_{x}f_{y}(1-s_{n})C\Psi\\
  (1-s^{2}f^{2}_{x}(1-s_{n}))C\Psi-s^{2}f_{x}f_{y}(1-s_{n})S\Psi\\
  s_{n}f_{x}S\Psi+s_{n}f_{y}C\Psi
  \end{bmatrix}
\end{equation}

As vectors $\hat{p}$,$\hat{e}$,$\hat{h}$ and $\hat{l}$ are co-planar as they hold common normal vector $\hat{n}$, let acute angle between $\hat{p}\perp\hat{e}$ and $\hat{h}\perp\hat{l}$ be $\xi$. Therefore,

\begin{equation}\label{xi}
\begin{array}{l}
     \cos(\xi)=\hat{l}.\hat{e}= sf_{y}sin\Psi-sf_{x}\cos\Psi,\\
     \sin(\xi)=|\hat{l}\times\hat{e}|= sf_{x}\sin\Psi+sf_{y}\cos\Psi
\end{array}
\end{equation}

Now, from the eq. (\ref{zetaddot3D}), the $\{\hat{p},\hat{e},\hat{n}\}$ frame with the angular accelerations $\ddot{\zeta}_{\hat{p}}$,$\ddot{\zeta}_{\hat{e}}$ and $\ddot{\zeta}_{\hat{n}}$ are produced by moments $\tau_{\hat{p}}$, $\tau_{\hat{e}}$ and $\tau_{\hat{n}}$, leads to:

\begin{equation}\label{penzddot}
\begin{bmatrix}
\ddot{\zeta}_{\hat{p}}\\
\ddot{\zeta}_{\hat{e}}\\
\ddot{\zeta}_{\hat{n}}
\end{bmatrix}=-(\Tilde{I}_{0}+\Tilde{I}_{R})^{-1}\begin{bmatrix}
\tau_{\hat{p}}\\
\tau_{\hat{e}}\\
\tau_{\hat{n}}
\end{bmatrix}-(\Tilde{I}_{0}+\Tilde{I}_{R})^{-1}\begin{bmatrix}
0\\
W_{||}R\\
0
\end{bmatrix}
\end{equation}

 Defining $\Tilde{I}$ as an overall inertia tensor matrix being $\Tilde{I} = \Tilde{I}_{0}+\Tilde{I}_{R}$. Where the $\Tilde{I}_{0}$ is the inertia tensor matrix calculated about the center of gravity (CG) and $\Tilde{I}_{R}$ is the inertia tensor matrix calculated about point of contact $P_{0}$ considering the robot as a point mass concentrated at CG. $W_{||}$ is the parallel component of the Weight vector $\Vec{W}$ along $\hat{p}$.

 \begin{equation}
     W_{||} = [0,0,-Mg]^{T}.\hat{p}=ss_{n}(f^{2}_{x}+f^{2}_{y})Mg
 \end{equation}

 Considering (\ref{xi}) the transform between $[\ddot{\zeta}_{\hat{p}},\ddot{\zeta}_{\hat{e}},\ddot{\zeta}_{\hat{n}}]^{T}$ and $[\ddot{\Theta},\ddot{\Phi},\ddot{\Psi}]^{T}$ is given by,

 \begin{equation}\label{Rxin}
  \begin{array}{l}
      [\ddot{\zeta}_{\hat{p}},\ddot{\zeta}_{\hat{e}},\ddot{\zeta}_{\hat{n}}]^{T} = R_{\xi\hat{n}}[\ddot{\Theta},\ddot{\Phi},\ddot{\Psi}]^{T},\\
      R_{\xi\hat{n}}=\begin{bmatrix}
      -sf_{x}S\Psi-sf_{y}C\Psi&-sf_{x}C\Psi+sf_{y}S\Psi&0\\
      -sf_{x}C\Psi+sf_{y}S\Psi&sf_{x}S\Psi+sf_{y}C\Psi&0\\
      0&0&1
      \end{bmatrix}
  \end{array}
 \end{equation}

Similarly, the relation between the moments $[\tau_{\hat{p}},\tau_{\hat{e}},\tau_{\hat{n}}]^{T}$ and $[\tau_{\hat{X}},\tau_{\hat{Y}},\tau_{\hat{Z}}]^{T}$ is given by,

\begin{equation}\label{TxyzTpen}
[\tau_{\hat{p}},\tau_{\hat{e}},\tau_{\hat{n}}]^{T}=R_{\hat{p}}R_{\alpha\beta\gamma}^{T}[\tau_{\hat{X}},\tau_{\hat{Y}},\tau_{\hat{Z}}]^{T},R_{\hat{p}}=[\hat{p},\hat{e},\hat{n}]
\end{equation}

Therefore, by combining the eq. (\ref{penzddot}),eq. (\ref{Rxin}) and eq. (\ref{TxyzTpen}) leads to the eq. (\ref{genangddot}).

\section{Inertia Tensor $\Tilde{I}$}
\label{'Inot'}

Consider a 3D object having mass $M$, undergoing general motion in the frame of reference $\mathscr{W}=\{\hat{x},\hat{y},\hat{z}\}$, the inertia tensor $\Tilde{I}$ of the object with respect to $\mathscr{W}$ is defined by:

\begin{equation}\label{inertiatensor}
\Tilde{I}=\begin{bmatrix}
\int(y^2+z^2).dM&-\int xy.dM&-\int zx.dM\\
-\int xy.dM&\int(z^2+x^2).dM&-\int yz.dM\\
-\int zx.dM&-\int yz.dM&\int(x^2+y^2).dM
\end{bmatrix}
\end{equation}

$\Tilde{I}$ depends on the geometry of the object. The matrix $\Tilde{I}$ is constructed with respect to specified origin and orthogonal axis. According to the parallel axis theorem, if inertia tensor $\Tilde{I}_{cm}$ is known about the CG of the object then $\Tilde{I}$ about any parallel axis is given by $\Tilde{I}=\Tilde{I}_{cm}+\Tilde{I}_{R}$ whereas $\Tilde{I}_{R}$ matrix is obtained by treating the object as a point mass at the CG.

\section*{Acknowledgment}

This work was supported by IvLabs, Visvesvaraya National Institute of Technology, Nagpur, India.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}