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

Sub-diffusive processes are present in several areas of natural sciences \cite{Silvestri2009,Weiss2007,Weiss2004,Scher1975}. The fundamental understanding of these phenomena is necessary due to the amounting experimental evidence of sub-diffusive behavior in complex physical, chemical and biological systems \cite{Caspi2000,Saxton2006,Golding2006}. A process is determined to be sub-diffusive when the {\it mean squared displacement} does not grow linearly in time, presenting a growth exponent smaller than one. This means that the process is significantly slower than normal diffusion. Also, by definition, sub-diffusion is a type of anomalous diffusion, a matter with great scientific interest.

The evolution of systems described by non-integrable Hamiltonians often exhibits regularity and chaos. The result is a mixed phase space containing chaotic seas, invariant tori, and Kolmogorov-Arnold-Moser (KAM) islands \cite{Lichtenberg1992}. For strongly chaotic systems,  the diffusion of a chaotic orbit through its phase space is considered to be normal, with no anomalous effects \cite{Ott1979,Zaslavsky2007}. For mixed-phase spaces, however, a chaotic orbit evolved from initial conditions near a KAM island may exhibits very complicated dynamical behavior, influenced by its neighborhood and other structures that interfere with the orbits transport. In this scenario, the diffusion is often anomalous \cite{Venegeroles2008,Diaz2019}.

Generally, the anomalous diffusion of orbits in the phase space of Hamiltonian systems is due to {\it stickiness} effect \cite{Zaslavsky2002,Altmann2013}. A chaotic orbit evolved from an initial condition set on a region surrounded by stability islands often experience {\it dynamical trappings}, spending a considerable amount of time in sticky regions around these islands and its cantori \cite{Zaslavsky1985,Balescu1997}. Applications of this trapping phenomenon can be found in many research areas as: fluid mechanics \cite{Solomon1993}, plasma physics \cite{delCastillo-Negrete2005}, celestial mechanics \cite{Contopoulos2010}, acoustics \cite{Altmann2009}, and biology \cite{Tel2005}. Additionally, studies on transport properties of chaotic orbits in Hamiltonian systems can be extended and applied to complex networks analysis \cite{Zou2016,Posadas-Castillo2014}.

In this work, we investigate the presence of sub-diffusive behavior in the Chirikov-Taylor Standard Map. First, we rely on the definitions of ergodicity and mixing to properly characterize the dynamic behavior of an ensemble of orbit evolved from an Unstable Period Orbit (UPO). The evolution started on the UPO guarantees stickiness influence in the chosen setup for the Standard Map. Based on the Mixing measure, we show that this dynamics can be described by a Continuous Time Random Walk (CTRW) model. In that sense, we chose to study the Fractional Diffusion Equation (FDE) to approach this problem. The solution of the FDE can be written in terms of the Mittag-Leffler function that depends explicitly on the anomalous diffusion exponent. Then, to be able to numerically calculate the anomalous diffusion exponent for the Standard Map, we choose a variant of the Ulam method. This method allows us to approximate the Perron-Frobenius operator for chaotic Hamiltonian maps. With that, it is possible to calculate the diffusion modes via an eigenvalue problem, and compare the computed eigenvalues to a fit provided by the Mittag-Leffler function, establishing a connection to the solution of the FDE. This comparison grants us an approximated value of the anomalous diffusion exponent that indicates a sub-diffusive behavior in the Standard Map, corroborating with other findings in the literature of stickiness and anomalous transport in Hamiltonian maps.

The paper is organized as follows. In Sec.\ \ref{sec:stdmap} we revise key aspects of the Standard Map. In Sec.\ \ref{sec:erg_mix} we briefly discuss ergodicity and import the definition of mixing to properly characterize the dynamical behavior of the selected setup for the Standard Map. In Sec.\ \ref{sec:anoma_diff} we review the CTRW model and discuss the analytical description of anomalous diffusion via the FDE. In Sec.\ \ref{sec:ulam_method} we present our selected numerical procedure to connect the dynamics in the Standard Map to anomalous diffusion, also referring to the Appendix where we further elucidate our numerical analysis. Finally, in Sec.\ \ref{sec:concluds} we conclude the chosen approach in this work.

\section{Standard Map}
\label{sec:stdmap}

The Standard Map \cite{Chirikov1969,Chirikov1979} can be used to describe the motion of a particle constrained to a movement on a ring while kicked periodically by an external field.

The map derives from the Hamiltonian
\begin{eqnarray}
H\left(q,p,t\right) & = & \frac{p^2}{2mr^2}+K\cos\left(q\right)\sum_{n=-\infty}^{\infty}\delta\left(\frac{t}{T}-n\right)~,
\label{eq:ham_stdmap}
\end{eqnarray}
where $\delta$ is the Dirac delta function, $q$ is the angular coordinate and $p$ is its conjugate momentum. It is worth remarking that we consider the particle's mass $m$, the ring's radius $r$, and the period of kicks $T$ all equal to one for simplicity.

Although the Hamiltonian in Eq.\ (\ref{eq:ham_stdmap}) is sufficient to analyse the dynamics of our interest, it is possible to define a sympletic non-linear discrete map $T_{SM}$ to investigate the dynamics via extensive and long numerical simulations, considering that iterating the sympletic map is much faster than solving the equations of motion.

The mapping $T_{SM}\left(p_{n},q_{n}\right)=\left(p_{n+1},q_{n+1}\right)$ gives the position and momentum for the $(n+1)^{th}$ iteration by the following equations
\begin{equation}
T_{SM}:\left\{\begin{array}{ll}
p_{n+1} = p_n +k\sin(q_n) \mod(2\pi)\\
q_{n+1} = q_n + p_{n+1} + \pi \mod(2\pi)\\
\end{array}
\right.,
\label{eq:stdmap}
\end{equation}
where the parameter $k$ controls the intensity of the non-linearity and the added term $+\pi$, on the equation for $q_{n+1}$, is to centralize the main island on the drawn phase spaces. It is also important to mind that this is an area-preserving map since the determinant of its Jacobian matrix is equal to unity. The phase space for the selected control parameter $k=1.46$ is drawn in Fig.\ \ref{fig:phase_space}.

\begin{figure}[h]
\includegraphics[scale=0.5]{Figs/phase_space_upo.png}
\caption{Phase space of the Standard Map for $k=1.46$. In blue the chosen initial conditions iterated until $10^5$ and in red the period $6$ Unstable Periodic Orbit (UPO).}
\label{fig:phase_space}
\end{figure}

One can notice in Fig.\ \ref{fig:phase_space} that the region between the main stability island and resonant period-6 satellite islands is denser than the outer chaotic sea region. This is a classic evidence of stickiness. In this case, orbits started from this region experience successive traps, filling this region more densely at first. Then, after a sufficient number of iterations, the orbits, once trapped, are free to access other regions of the chaotic sea.

In Fig.\ \ref{fig:phase_space} we are also drawing in red a period-6 UPO, associated with the period-6 satellite stable islands. The position of this UPO was calculated numerically with a floating-point double-precision via the method outlined in \cite{Ciro2018}. It is important to properly determine the position of this selected UPO because it will set our ensemble of initial conditions for the upcoming analysis.

\section{Ergodicity and Mixing}
\label{sec:erg_mix}
In physics, ergodicity is commonly defined by the equivalence of the ensemble and time averages of a relevant observable $y$, provided such averages formally exists:

\begin{equation}
    \langle y\rangle_{ens}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_0^{T}y(t^{\prime})dt^{\prime}~.
    \label{eq:ergodicity}
\end{equation}Since the right-hand side cannot depend on time, the process has to be stationary (so that the left-hand side is
time-independent) and possess an invariant probability measure, over which the ensemble mean in the left-hand side
is calculated.

Experimentally\footnote{In our case, ``experimentally'' means in numerical simulations written in JULIA language.} this means that if Eq. (\ref{eq:ergodicity}) holds, averaging many single-time measurements over repeated experiments on the ensemble of similar systems, yields the same results as a time average over a single, but prolonged experimental run \cite{Meroz2015}. In experimental data, it is necessary to use a discrete-time variant of Eq. (\ref{eq:ergodicity}) given by

\begin{equation}
    \langle y\rangle_{ens}=\lim_{N\rightarrow\infty}\frac{1}{N}\sum^{N-1}_{n=0}y_n~.
\end{equation}

Mixing, however, is the asymptotic independence of $y_n$ and $y_0$, as the step number $n$ goes to infinity. Here $y_n$ is our relevant observable and it behaves as random variable. Mixing is a stronger property than ergodicity \cite{Lasota1994} because all mixing processes are ergodic, but not vice versa. Proving mixing is enough to show ergodicity, and in some cases this may be an easier task.

In \cite{Garbaczewski1995,Janicki1993}, a 1D dynamical functional is defined by

\begin{equation}
    D_n=\langle e^{i(y_n-y_0)} \rangle~.
    \label{eq:dyn_func}
\end{equation}With that, the dynamical functional test \cite{Magdziarz2011} is based on the fact that $E_n = D_n - |\langle e^{i y_0} \rangle|^2$ vanishes when the random variables $y_n$ and $y_0$ are independent. Thus, a stationary process $y_n$ is said to mix if and only if

\begin{equation}
    \lim_{n\rightarrow \infty}E_n=0~.
    \label{eq:mixi}
\end{equation}

To show that the process is mixing, ensemble averages need to be calculated. Moreover, the process is ergodic if and only if

\begin{equation}
    \lim_{n\rightarrow \infty}\frac{1}{n}\sum_{j=0}^{n-1}E_j=0~.
    \label{eq:erg_mix}
\end{equation}

Note that the tests are quite restrictive with respect to the amount of data, which has to be sufficiently large.

\begin{figure}[h]
\includegraphics[scale=0.5]{Figs/erg_mix_pq.png}
\caption{Mixing measure for the $p$ and $q$ series, considering an evolution started at the period-6 UPO for the selected value of the control parameter $k=1.46$ of the Standard Map.}
\label{fig:erg_mix_series}
\end{figure}

For the coordinates $q$ and $p$ of the Standard Map, we make use of the mixing definition considering $y_n=p_{n+\Delta n}-p_n$ and $y_n=q_{n+\Delta n}-q_n$, as the increment process for each coordinate. We evolved an ensemble of $10^3$ initial conditions in a small square of size $10^{-4}$, centered at the position of the period-6 UPO calculated before and, with that we analyse the time-series of $q$ and $p$. We show in Fig.\ \ref{fig:erg_mix_series} the asymptotic behavior of $E_n$ for $p$ and $q$ series.

From Fig.\ \ref{fig:erg_mix_series} it is clear that either $p$ or $q$ does not present ergodic behavior, since Eqs.\ (\ref{eq:mixi},~\ref{eq:erg_mix}) do not hold. This is expected because the selected setup for the Standard Map has a mixed phase space, evidenced in Fig.\ \ref{fig:phase_space}. A mixed phase space contains regions of stability along with a chaotic sea, however, due to this mix, other relevant structures that influence the transport of orbits are also present, this is the case of invariant manifolds associated with the dynamics. It is also important to note that we chose the ensemble centered at the UPO, guaranteeing a dynamics influenced not only by the stickiness but also by the unstable and stable branches of the manifold associated with this UPO. An ensemble evolved from other regions of the chaotic sea may present a dynamic behavior that agrees with Eqs.\ (\ref{eq:mixi},~\ref{eq:erg_mix}).

Furthermore, the results of the mixing measure for the dynamics of the Standard Map, shown in Fig.\ \ref{fig:erg_mix_series}, allow a different interpretation for this dynamics. According to the decision tree proposed in \cite{Meroz2015}, since the trajectories of our interest exhibit a non-ergodic dynamical behavior, they can be modeled by the Continuous Time Random Walk (CTRW) description. We explain the CTRW model in detail and connect it to a theoretical description of anomalous diffusion in the next section.

\section{Continuous Time Random Walk and Anomalous Diffusion}
\label{sec:anoma_diff}

The CTRW model was introduced by Scher and Montroll when they investigated the anomalous transport properties of charge carriers in amorphous materials \cite{Scher1975}. The model is given by successive traps of a random walk, in which each step is characterized by two independent stochastic processes, the waiting time and the displacement in space. The combination of these two factors yields an anomalous diffusion process. It is important to note that, because of the basic nature of CTRW, it is a well-suited candidate to describe stickiness phenomena in Hamiltonian systems.

Since CTRW is an example of anomalous diffusion, we studied the Fractional Diffusion Equation (FDE), one of the main theoretical descriptions of anomalous diffusion. Further discussions and examples of CTRW simulations are presented in the Appendix. The FDE is given by

\begin{eqnarray}
\frac{\partial^{\alpha}\rho\left(x,t\right)}{\partial t^{\alpha}} & = & D\frac{\partial^{2}\rho\left(x,t\right)}{\partial x^{2}}~,
\label{eq:FDE}
\end{eqnarray}note the presence of the operator $\frac{\partial^{\alpha}}{\partial t^{\alpha}}$, (Caputo fractional derivative with $0<\alpha\leq1$) \cite{Haubold2011}. The parameter $\alpha$ represents the anomalous diffusion exponent, and $D$ the
diffusion coefficient.

In a closed interval $0\leq x\leq 2\pi$ with periodic boundary conditions, the solution of Eq. (\ref{eq:FDE}) is given by the following superposition \cite{Haubold2011}

\begin{equation}
   \text{\small$
   \rho\left(x,t\right) =  a_0+\sum_{m=1}^{\infty}E_{\alpha}\left(-D m^{2}t^{\alpha}\right)\left(a_{m}\cos\left(mx\right)+b_{m}\sin\left(mx\right)\right)~,$}
   \label{eq:supsolFDE}
\end{equation}where the first term in the right-hand side can be set to $a_0=\frac{1}{2\pi}$ due to normalization of
$\rho\left(x,t\right)$. The terms $a_{m}$ and $b_{m}$ are related to the initial condition $\rho\left(x,0\right)$ and $E_{\alpha}\left(z\right)$ is the Mittag-Leffler function defined by

\begin{eqnarray}
E_{\alpha}\left(z\right) & = & \sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma\left(1+\alpha k\right)}~.
\label{eq:MittagL}
\end{eqnarray}

Here $\alpha, z \in \mathbb{C}$, $\Re\left(\alpha\right)>0$, $\mathbb{C}$ being the set of complex numbers and  $\Gamma\left(z\right)$ is the Gamma function \cite{Haubold2011}.

Regarding the solution presented in Eq.\ (\ref{eq:supsolFDE}) and the Mittag-Leffler function defined in Eq.\ (\ref{eq:MittagL}) the following observations can be made:

\begin{enumerate}
\item The density $\rho(x,t)$ approaches to a limit $\rho\sim\frac{1}{2\pi}$;

\item Besides of the first term $a_0$ of Eq. (\ref{eq:supsolFDE}), all other terms of the summation have zero integral in space, e.g. $\int_{0}^{2\pi}\cos\left(mx\right)dx=0$, etc;

\item The Mittag-Leffler function equals to the exponential when $\alpha=1$;

\item For values of $0<\alpha<1$ the decay given by $E_{\alpha}\left(-t^{\alpha}\right)$ is slower than exponential but the limit $\lim_{t\rightarrow\infty}E_{\alpha}\left(-t^{\alpha}\right)=0$ still holds.
\end{enumerate}

Furthermore, one can define

\begin{eqnarray}
\text{\small$
\hat{M} \left(x,y,t\right) = \frac{1}{2\pi}+\frac{1}{\pi}\sum_{m=1}^{\infty}E_{\alpha}\left(-D m^{2}t^{\alpha}\right)\left[\cos\left(my\right)\cos\left(mx\right)+\sin\left(my\right)\sin\left(mx\right)\right]~,$}
\label{eq:Op_FDE}
\end{eqnarray} and rewrite Eq. (\ref{eq:supsolFDE}) as following

\begin{eqnarray}
\rho\left(x,t\right) & = & \int_0^{2\pi} dy ~\hat{M} \left(x,y,t\right) \rho\left(y,0\right)~.
\label{eq:TransOp_flow}
\end{eqnarray}

Then, Eq.\ (\ref{eq:TransOp_flow}) defines a Perron-Frobenius-like transfer operator by virtue of the properties of the solutions $\rho\left(x,t\right)$ in Eq. (\ref{eq:FDE}). In one hand, the eigenfunctions of the operator, given by $\left[1,~\cos\left(mx\right),~\sin\left(mx\right)\right]$, $m\in\mathbb{Z}^{+}$ are time-independent. On the other, the eigenvalues, given by $\left[1,~E_{\alpha}\left(-Dm^{2}t^{\alpha}\right)\right]$, carry the time dependence.

One can observe that changing the right-hand side of Eq. (\ref{eq:FDE}) for other linear position-dependent operator, only changes the eigenfunctions and not the eigenvalues. Also, the presence of $m$ in the argument of the Mittag-Leffler function  makes the influence of the eigenfunctions with higher $m$ go to zero fast. Then, it is fair to assume that, in the asymptotic time $t\rightarrow\infty$, the transient behavior is governed by the first eigenfunction.

A discrete scenario of Eq.\ (\ref{eq:TransOp_flow}) can be used to analyze the dynamics governed by mappings. The discrete version of Eq.\ (\ref{eq:TransOp_flow}) is given by

\begin{eqnarray}
\rho\left(i,n\right) & = & \sum_{j} ~\hat{S} \left(i,j,n\right) \rho\left(j,0\right)~,
\label{eq:TransOp_disc}
\end{eqnarray} where $i,j$ are {\it space-like} indexes and $n$ a {\it time-like}. Assuming that $[\lambda_l,\psi_l]$ are the eigenpairs associated to $\hat{S}\left(i,j,n\right)$ and, since this matrix is a Perron-Frobenius-like operator, the following remarks are noted:

\begin{enumerate}
\item The index $l$ numerates the eigenpairs $[\lambda_l,\psi_l]$;
\item All eigenvalues of $\hat{S}\left(i,j,n\right)$  have the constrain $\left|\lambda_{l}\left(n\right)\right| \leq 1$;
\item There is at least one eigenvector $\psi_{0}\left(i\right)$ with eigenvalue $\lambda_{0}=1$ and $\sum_{i}\psi_{0}\left(i\right)=1$;
\item Eigenvectors $\psi_{l}$ associated to $|\lambda_l|<1$, have  $\sum_{i}\psi_{l}\left(i\right)=0$. Otherwise the probability would not be conserved.
\end{enumerate}

For a Markovian process ${\footnotesize\hat{S}\left(i,j,n\right)=\left(\hat{S}\left(i,j,1\right)\right)^{n}}$, and the eigenvalues of $\hat{S}\left(i,j,n\right)$ are given by $\lambda\left(n\right)=\lambda^{n}$, where $\lambda$ represents the eigenvalues of $\hat{S}\left(i,j,1\right)$. In the case that there is just one eigenvalue with $\lambda_{0}=1$, the time dependence of the probability distribution is given by

\begin{eqnarray}
\rho\left(i,n\right) & = & \psi_{0}\left(i\right)+\sum_{l}\exp\left(n\log\left(\lambda_{l}\right)\right)a_{l}\left(0\right)\psi_{l}\left(i\right)~,
\label{eq:normaldiffsol}
\end{eqnarray} one can note that, in this case, the solution corresponds to the superposition solution shown in Eq.\ (\ref{eq:supsolFDE}) with $\alpha=1$.

For a non-Markovian process in general ${\footnotesize\hat{S}\left(i,j,n\right)\neq\left(\hat{S}\left(i,j,1\right)\right)^{n}}$ and the solution given by Eq. (\ref{eq:normaldiffsol}) is not longer valid. Nevertheless, the following behavior shown in Eq. (\ref{eq:eigenvalue_Mittag}) can be valid if the diffusion process is believed to come from a CTRW-like model
\begin{eqnarray}
\lambda_l\left(n\right) & = & E_{\alpha_l}\left(-d_l n^{\alpha_l}\right)~,
\label{eq:eigenvalue_Mittag}
\end{eqnarray}where we define $\alpha_l$ and $d_l$ as, respectively, the anomalous diffusion exponent and coefficient related to the $l^{th}$eigenpair, see the Appendix for further discussion and examples.

Our approach in this work is focused in the eigenpair with the highest eigenvalue smaller than one $[\lambda_1,\psi_1]$, since, as stated before, it is the one that governs the transient dynamics of $\rho$ in the asymptotic limit  $n\rightarrow\infty$.

Finally, is important to note from Eq. (\ref{eq:eigenvalue_Mittag}), that it is possible to assess the anomalous diffusion exponent $\alpha_1$, associated to the assumed most relevant eigenpairs $[\lambda_1,\psi_1]$. Hence, the problem now is to calculate the Perron-Frobenius-like operator for the selected setup of the Standard Map to characterize the most relevant dynamic behavior via $\alpha_1$. In the next section we discuss a numerical approach to this problem.

\section{Numerical approach - Ulam method}
\label{sec:ulam_method}

\begin{figure}[b]
\includegraphics[scale=0.28]{Figs/diff_mode_.png}
\caption{First diffusion mode, given by the first eigenvector $\psi_1$ of $\hat{S}\left(n^{\prime}\right)$, on the phase space of the standard map for the selected control parameter $k=1.46$. We considered the space divided by a grid of $512\times512$ cells.}
\label{fig:diff_mode}
\end{figure}

To be able to simulate the discrete version of the transfer operator equation discussed before, considering the dynamics of the selected setup for the Standard Map, we make use of a variant of the Ulam method proposed by Frahm and Shepelyansky in \cite{Frahm2010,Frahm2013}.

The chosen numerical method is used to approximate the Perron-Frobenius operator for chaotic Hamiltonian maps. It is based on following the evolution of one chaotic orbit through a discretized phase space. First, we need to divide the phase space into a regular grid and, at each iteration, ask, in which cell $j$ the orbit is, and at which cell $i$ it jumps in the next iteration. After many iterations, we calculated how many times an orbit in cell $j$ jumped to cell $i$. Then, it is defined the matrix elements $S_{ij}$. Considering a normalization condition $\sum_{i}S_{ij}=1$, that give us the probability to jump from $j$ to $i$ after one iteration ($n^{\prime}=1$). Additionally, if we consider the cell $i$ that we arrive from cell $j$ after $n^{\prime}$ iterations we have the matrix element $S_{n^{\prime},ij}$. Due to the normalization condition, the matrices $\hat{S}\left(n^{\prime}\right)$ with elements $S_{n^{\prime},ij}$, are Perron-Frobenius-like operators, which have at least one invariant probability vector. It is important to note that we differentiate $n$ from $n^\prime$ as, respectively, the number of iterations for the orbit's evolution and the iterations skipped to calculate the Perron-Frobenius-like operator.

The eigenvectors, in this case, are called diffusion modes and are defined by the eigenvalue problem

\begin{eqnarray}
\sum_{j}S_{ij}\psi_{l}\left(j\right) & = & \lambda_{l}\psi_{l}\left(i\right)~,
\end{eqnarray}
where $\lambda_{l}$ is the eigenvalue and $\psi_{l}\left(i\right)$ is the value of the eigenvector $\psi_{l}$ at the cell $i$. Note that $\psi_{l}$ may have different numerical values for different cells in the correspondent region of the discrete space. We call the first diffusion mode the $\psi_{l}$ that has the largest $\left\Vert \lambda_{l}\right\Vert \neq1$.

In discrete spaces is often necessary to set a fine grid that guarantees enough resolution to distinguish small structures. However, a fine grid combined with sufficient $n^{\prime}$ iterations to calculate the matrices, demand powerful computational resources. To surpass this issue, we make use of the Arnoldi method described in \cite{Frahm2010,Frahm2013}. With that, we are able to approximate the eigenpairs $\left[\lambda_l,\psi_l \right]$.

In Fig.\ \ref{fig:diff_mode} we show the numerical result for the first diffusion mode $\psi_1$ through the phase space of the Standard Map, considering only one trajectory started from the selected period-6 UPO, but evolved until $10^{10}$ iterations. For this first result we considered $n^{\prime}=1$, since $\psi_1$ does not change enough to modify this analysis. This matter will be addressed further in this section.

\begin{figure}[b]
\includegraphics[scale=0.42]{Figs/fit_stdmap_.png}
\caption{Eigenvalue $\lambda_{1}$ of matrix $\hat{S}\left(n^{\prime}=1\right)$ after skipping $n^{\prime}$ iterations for the Standard map with $k=1.46$. The analytical fitting was provided by a numerical representation for the Mittag-Leffler function.}
\label{fig:fit_std_map}
\end{figure}

It is clear from the color map in Fig.\ \ref{fig:diff_mode} that the phase space is divided in two distinct regions. This division is properly characterized by the change of sings of the first diffusion mode $\psi_1$. Negative values of $\psi_1$ are attained only at the region around the main island and the other resonant ones. Yet, positive values of $\psi_1$ depicts the chaotic sea.
The regions with $\psi_1=0$ are also interesting, but will be investigated in another study.

One can notice that the region of negative $\psi_1$ is close related to the region of stickiness influence, as depicted early in Fig.\ \ref{fig:phase_space}. This establishes a first connection between the dynamics affected by the stickiness and its diffusion mode. To further investigate this special region, we analyse the behavior of the first eigenvalue $\lambda_1$, associated with the first diffusion mode, as function of the first nine iterations skipped to calculate the Perron-Frobenius-like operator. Moreover, based on the analytical description made in last section, we propose fitting the numerical values of $\lambda_1$ with a numerical representation \cite{Gorenflo2002,Valerio2014} for the Mittag-Leffler function at Eq.\ (\ref{eq:MittagL}). This result is shown in Fig.\ \ref{fig:fit_std_map}.

The connection between the computed eigenvalues $\lambda_1=\lambda_1(n^{\prime})$ and the fit provided by the Mittag-Leffler function is further addressed and tested on the Appendix. There we simulate the CTRW model for different values of the anomalous diffusion coefficient $\alpha$ and show the reliability of this procedure.

\begin{figure}[b]
\includegraphics[scale=0.42]{Figs/error_fit_std.png}
\caption{Dependence of the error function defined in Eq.\ (\ref{eq:error_fit}) with the number of iterations skipped, i. e. the dot product between the diffusion mode after one iteration and the diffusion mode after $n^{\prime}$ iterations. The dashed line at $\err(\psi,1,n^{\prime})=1$ indicates the ideal value.}
\label{fig:error_fit_std_map}
\end{figure}

The result shown in Fig.\ \ref{fig:fit_std_map} provides a fitted value of $\alpha_1=0.250\pm0.008$ that implies a dependence slower than exponential for $\lambda_{1}(n^{\prime})$. This corroborates our initial assumption that the dynamic behavior of this setup for the Standard Map is indeed related to anomalous diffusion, specifically a sub-diffusive process with $\alpha\approx0.25$.

It is also important to mention that, when studying the dependence of $\lambda_{1}$ with the matrix iteration $n^{\prime}$, we must check if we are still considering the same diffusion mode, or, in each iteration, it is changing due numerical errors. To control this deviance from unity, we apply the normalization conditions $\sum_{i}\psi_{1}\left(i,t=1\right)\psi_{1}\left(i,t=1\right)=1$, $\sum_{i}\psi_{1}\left(i,t=n^{\prime}\right)\psi_{1}\left(i,t=n^{\prime}\right)=1$
and compute the value of the projection
\begin{equation}
P(\psi,1,n^{\prime})=\left\Vert\psi_{1}\left(t=1\right)\cdot\psi_{1}\left(t=n^{\prime}\right)\right\Vert=\left\Vert \sum_{i}\psi_{1}\left(i,t=1\right)\psi_{1}\left(i,t=n^{\prime}\right)\right\Vert~,
\label{eq:error_fit}
\end{equation}
i.e., the dot product between the diffusion mode after one iteration and the diffusion mode after $n^{\prime}$ iterations.

In Fig.\ \ref{fig:error_fit_std_map} we show the computed projection $P(\psi,1,n^{\prime})$ for the first nine iterations skipped to calculate the Perron-Frobenius-like operator. Since the error does not diverge from its ideal value of unity, considering only the first decade, it is fair to assume that our analysis of $\alpha\approx0.25$ still holds. However, for larger $n^{\prime}$, we are not dealing with the same diffusion mode, and the calculation of $\alpha$ becomes meaningless.

It is important to note that with the described numerical procedure is possible to calculate the most relevant diffusion mode $\psi_1$, and associated with it, the anomalous diffusion exponent $\alpha_1\approx\alpha$ that characterize the predominant dynamical behavior for the time evolution of interest. In the case of the selected setup for the Standard Map, the anomalous sub-diffusive behavior is confirmed by $\alpha\approx0.25$. Furthermore, one can adapt the procedure to investigate other dynamical systems, particularly to study sub-diffusive transport in complex networks.

\section{Conclusions}
\label{sec:concluds}

In this work, we present a procedure, based on different descriptions found in the literature of anomalous diffusion, to show that the evolution of trajectories under stickiness influence in the Chirikov-Taylor Standard Map can indeed be described as a sub-diffusive process.

Initially, we showed that for the chosen control parameter $k=1.46$, the Standard Map presents a non-ergodic dynamical behavior, as expected for a mixed phase space. Then, we assumed that the trajectories in this setup can be described by the Continuous Time Random Walk model. Since CTRW is a classic example of anomalous diffusion, we studied some consequences of the Fractional Diffusion Equation and how to connect it to a numerical method that approximates the Perron-Frobenius operator for the map. With that, we established a relation between the eigenvalues of the Perron-Frobenius operator and the solution of FDE, providing an approximated value for the anomalous diffusion exponent $\alpha$ for the selected setup of the Standard Map. The value of $\alpha\approx0.25$ shows that the evolution of trajectories in this scenario is indeed associated with anomalous diffusion, distinctively a sub-diffusive behavior. Moreover, we connected the stickiness in the Standard Map with a robust framework to describe this effect as a sub-diffusive anomalous transport.

It is important to note that the procedure described here is readily applicable to other Hamiltonian systems with a variety of applications and also suitable to investigate transport properties on complex networks.

\section*{Acknowledgments}
This study was financed in part by the Coordena\c c\~ao de Aperfei\c coamento de Pessoal de N\'ivel Superior - Brazil
(CAPES) - Finance Code 001, S\~ao Paulo Research
Foundation (FAPESP) - Brazil, under Grants No. 2018/03211-6 and 2018/03000-5, and IRTG 1740 financed by Deutsche Forschungsgemeinschaf (DFG).

\section*{Appendix: Stochastic representation and RW/CTRW Simulations}
\label{app:A}

This Appendix is devoted to verify and support our analysis made for the selected setup of the Standard Map, based on the connection between the eigenvalue of the Perron-Frobenius transfer operator and the Mittag-Leffler function.

We address here how to simulate a Random Walk (RW) and a Continuous Time Random Walk (CTRW) using the It\^o stochastic representation \cite{Magdziarz2007}.

A classical CTRW is a process subordinated to a simple RW, in which steps (jumps) follow at random instants of time. The waiting time $t_i$ for the next step follows the probability distribution with the known probability density $\psi(t_i)$. The number of steps performed up to time $t$, $\nu(t)$ is the {\it operational time}, or the subordinator of the corresponding subordination scheme. The {\it clock time} $t$ of the $\nu$-th step is then the following sum
\begin{equation}
   t(\nu) = \sum_{i=1}^\nu t_i~,
   \label{eq:time_sum}
\end{equation}which, for fat-tailed $\psi(t_i)\sim t^{-1-\alpha}$
with $0< \alpha< 1$ tends in distribution to a one-sided Levy law. The value of $\nu$ as a function of $t$ is then $\nu(t) = \mathrm{inf}\{ \nu, t(\nu)>t \}$ \cite{Sokolov2005}. This gives us the prescription for stochastic simulations of the CTRW.

For long times the variable $\nu$ can be taken continuous; the corresponding continuous limit for the Probability Density Function $\rho(x,t)$ is
then given by the Fractional Diffusion Equation showed in Eq.\ (\ref{eq:FDE}).

In the displacement $X(t)$ then follows from a couple of stochastic differential equations \cite{Fogedby1994}
\begin{eqnarray}
dX(\nu) &=& \sqrt{2D}dB(nu), \\
dt(\nu) &=& \tau d L_\alpha(\nu)
\label{eq:displace}
\end{eqnarray}with $B(\nu)$ being a standard Brownian motoin, $L_\alpha(\nu)$ a one-sided L\'evy flight, and $D$ and $\tau$
the appropriate constants. From this representation Eq.\ (\ref{eq:FDE}) follows as well.

With that, we followed the numerical method proposed by Magdziarz and Weron in \cite{Magdziarz2011,Magdziarz2007} to simulate the dynamics of a RW process, setting $\alpha=1$ and CTRW processes, setting $\alpha=0.6$ (relatively not so long trap-time) and $\alpha=0.3$ (relatively long trap-time). In Fig.\ \ref{fig:ctrw_examples} we show examples of simulated trajectories of these three different dynamic scenarios.

\begin{figure}[h]
\includegraphics[scale=0.45]{Figs/ctrw_examples.png}
\caption{Examples of trajectories for a RW process ($\alpha=1)$ in black. Also, for two distinct CRTW processes considering a relatively not so long trap-time ($\alpha=0.6$) in red, and a long trap-time ($\alpha=0.3$) in green. The time-axis was chosen to be linear to properly illustrate the different trap-times.}
\label{fig:ctrw_examples}
\end{figure}

The discussed Ulam method, shown in Sec.\ \ref{sec:ulam_method}, used for approximate the Perron-Frobenius-like operator to approach the anomalous diffusion exponent $\alpha_1$ for the selected setup of the Standard Map, can also be applied to an ensemble of particles whose movement is governed by RW. However, it is not an appropriate method to analyze cases where $\alpha\neq1$. In that case of anomalous diffusion, it is necessary to apply a kernel density estimator to have a good approach to the probability density from a histogram \cite{Magdziarz2007}. To avoid the expensive use of a kernel density estimator, we choose to use a different method called the Extended Dynamic Mode Decomposition (EDMD) \cite{Willams2015,Klus2016}.

\begin{figure}[b]
\includegraphics[scale=0.45]{Figs/ctrw_fiting.png}
\caption{Behavior of the eigenvalue $\lambda_{1}$ as function of the log of the number of iterations $n^{\prime}$ to calculate the correspondingly matrices, considering the dynamics of a RW and the two CTRW examples. The points mark the simulated dynamics in each case and the dashed line is the fitting provided by a numerical representation of the Mittag-Leffler function.}
\label{fig:ctrw_fitting}
\end{figure}

The EDMD in \cite{Klus2016} is method to approximate the Perron-Frobenius operator, using the fact this operator and the Koopman operator\footnote{The Koopman operator maps functions of state space to functions of state space.} are adjoint to each other. The EDMD method focus on a dictionary of observables $D=\left[\phi_1(x),\phi_2(x),\ldots,\phi_k(x)\right]$, functions of state space, and how they change along the trajectory. With this information, and the definition of the Koopman operator is possible to approximate the eigenfunctions and eigenvalues both of the Koopman operator as of the Perron-Frobenius operator. Much of the method relies in an educated guess for the dictionary of observables, that must be rich enough to approximate the Koopman operator eigenfunctions, see \cite{Willams2015,Klus2016} for further details.

We consider the region $0\leq x\leq2\pi$ with periodic boundary conditions and choose the observables
\begin{equation}
   \text{\small$D= \left[1,\cos\left(x\right),\sin\left(x\right),\cos\left(2x\right),\sin\left(2x\right),\ldots,\cos\left(5x\right),\sin\left(5x\right)\right]$}~,
   \label{eq:observ}
\end{equation}as the dictionary for the EDMD method.

The behavior of the eigenvalue $\lambda_{1}$ as function of the skipped iterations $n^{\prime}$ is shown in Fig.\ \ref{fig:ctrw_fitting} for the same simulations depicted in Fig.\ \ref{fig:ctrw_examples}. Considering its relation to the solution of the FDE, as explained in Sec.\ \ref{sec:anoma_diff}, we selected the Mittag-Leffler function, given by

\begin{eqnarray}
E_{\alpha}\left(z\right) & = & \sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma\left(1+\alpha k\right)}~,
\label{eq:MittagL_}
\end{eqnarray}to fit these three different dynamic scenarios for three different values of $\alpha$.

The Mittag-Leffler is the suitable function because of the assumed behavior

\begin{eqnarray}
\lambda_l\left(n\right) & = & E_{\alpha_l}\left(-d_l n^{\alpha_l}\right)~,
\label{eq:eigenvalue_Mittag_}
\end{eqnarray}where it was defined $\alpha_l$ and $d_l$ as, respectively, the anomalous diffusion exponent and coefficient related to the $l^{th}$eigenpair of the Perron-Frobenius-like operator.

Is important to mind that in these particular cases, the anomalous diffusion exponent $\alpha_l$ related to the $l^{th}$eigenpair, is the same as $\alpha_1$. Then, the anomalous diffusion exponent $\alpha$, that characterize the anomalous dynamics, is $\alpha\approx\alpha_1=\alpha_l$.

Finally, one can observe that the values of the fitted functions are close to the actual values set for its simulations. Specially, for the CTRW model, the given values of $\alpha_1=0.637\pm0.009$ and $\alpha_1=0.329\pm0.004$ imply a dependence slower than exponential, as expected from the theory.


\section{INTRODUCTION \label{sec:intro}}
Research in the field of layered materials has grown quickly since the development of the high-quality sample yielding scotch-tape exfoliation technique.~\cite{Geim2013,Toth2019} For instance, the quantum spin Hall effect at up to 100 K was observed in monolayer tungsten ditelluride (WTe$_2$)~\cite{Wu76}. Another important discovery was the observation of a correlated insulator state by tuning the twist degree of freedom in bilayer structures of graphene, due to the presence of flat bands near zero Fermi energy at a twist angle of about $1.1^\circ$, which upon electrostatic doping yields superconducting states with a critical temperature of up to 1.7 K.~\cite{Cao2018corr,Cao2018unc} Research in the field of point defect qubit candidates has also grown rapidly, particularly since the detection of single negatively charged nitrogen vacancy (N$V^-$) color centers in diamond.~\cite{Gruber2012} These qubit candidates consist of defects in the crystal structure involving substitutional or interstitial atoms and/or vacancies, and act as single photon sources as well as sources of electronic spin. Desirable characteristics for such qubits include indistinguishability of emitted photons, negligible spectral diffusion and long spin coherence time. State-of-the-art results include second-long coherence time of N$V$ centers,~\cite{Gill2013} km-scale entanglement of two N$V$s,~\cite{Hensen2015}, and the discovery of spectrally stable germanium vacancy (Ge$V$) and silicon vacancy (Si$V$) point defects.~\cite{Siyushev2017,Becker2018} The properties of spin qubits can be tuned by changing the crystal structure and the nature of the defect.
Engineering spin qubits with properties that are desirable for applications in sensing, quantum computing and quantum communications are at the forefront of research interests.
Nevertheless, identifying the right combination of crystal structure and defect is enormously challenging due to the combinatorially large number of potential candidates that need to be explored using experiments or first-principles calculations.

The layered material hexagonal boron nitride ($h$-BN) is of particular interest among layered materials as it has a wide bandgap,~\cite{Elias2019} enabling its ability to encapsulate other layered materials and to host a variety of point defects giving rise to optical transitions lying within the range of its bandgap.~\cite{Toth2019,Tran2,IvadyhBN,Mehdi2018,Gottscholl2020} The issue with $h$-BN as a host for a qubit candidate employing electronic spins is that the atomic nuclei of boron and nitrogen have spins as well, which can cause spin-decoherence of the electronic spin state. Such an argument invites consideration of SiS$_2$ as a host, which theory has predicted to exist in layered form between 2.8 GPa and 3.5 GPa for the space group P$2_1$/c.~\cite{Plasienka2016} The SiS$_2$ host would be diamond-like or silicon carbide (SiC)-like in that the atoms that constitute it have a low natural abundance of isotopes with nuclear spin, which can be further suppressed via isotope purification in growth, diamond and SiC being systems that have been extensively studied as hosts for qubit candidates.~\cite{Kuate8,Becker2018,DOHERTY20131,Siyushev2017,Gali2010,SangYun2013,Weber8513,Soltamov2012,Koehl2011,Kraus2014,Wahl,Wang_2019,Dong,Green_2017,Green_2019,Kuate,DEFO2019,sohn2018cont} The point defect we will investigate in SiS$_2$ has the further advantage of exhibiting inversion symmetry, which eliminates the issue of an electric dipole moment making it susceptible to external noise and local fields and thus causing broadening of transitions.

Given the growing number of experimental investigations, the need to rapidly identify promising point defect qubit candidates is being recognized. Indeed, recent work~\cite{davidsson2021adaq} has looked into automating the process of characterizing point defects through a code named ADAQ, to aid in the more rapid identification of defect signatures from experiments, and a push toward high-throughput point defect calculations was also present in the earlier PyCDT code.~\cite{BROBERG2018165} In identifying point defects in experiments, knowledge of the zero-phonon line (ZPL) transition is essential. We present in this paper two methods that can be incorporated into high-throughput searches for point defect qubit candidates. The first is a quick new method for estimating the error associated with computing the energy of the zero-phonon line (ZPL) transition of a point defect using Janak's theorem and the second is a new method for rapid convergence of excited state self-consistent field (SCF) calculations. These new tools are an important step in accelerating the high-throughput discovery of materials that can support spin qubits for applications in sensing, quantum communication and quantum computing.

Janak's theorem states that,~\cite{Janak1978}
\begin{equation}
\frac{\partial E}{\partial n_i} = \epsilon_i,
\end{equation}
where $E$ is the total energy, $n_i$ is the orbital occupation of the $i^{\rm th}$ orbital, and $\epsilon_i$ is the corresponding eigenvalue of the orbital. This theorem is similar to Koopmans' theorem for Hartree-Fock theory~\cite{KOOPMANS1934} in that it relates energy differences under a change in the number of electrons to orbital energies. The theorem is in fact a density-functional theory (DFT) version of a theorem originally proved by Slater for his $X-\alpha$ method, which was introduced in an early attempt to account for both exchange and correlation in electronic structure calculations.~\cite{Slater1970,Kaxiras} As Janak's theorem does not imply that $\Delta E = \epsilon_i \Delta n_i$ for finite $\Delta n_i$, an estimate of the error associated with using the theorem to determine excitation energies where orbital occupations change by integer amounts is necessary.
Applying the theorem and the method of error estimation is inherently faster than computing the energy of the excited state for a given point defect, as the error calculation amounts to terminating the excited state calculation before full convergence. We also show that excited state calculations for systems with small strain converge faster if the charge density is initialized by appropriately mixing the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) from the ground state calculation, leaving the contribution to the charge density from the remaining orbitals unchanged, as compared to initialization of the charge density from a superposition of atomic charges.

In this work, we first outline the computational methods in Section \ref{sec:disc}. In Section \ref{sec:res} we include results related to (i) ZPL estimations for silicon monovacancies in $4H$-SiC (see Section \ref{sec:VSi}) and (ii) the singly negatively charged calcium vacancy in SiS$_2$ (see Section \ref{sec:CaV}, which also includes stability calculations) and a discussion of these results, ending with a summary of our conclusions.

\section{THEORETICAL FORMULATION AND APPROACH \label{sec:disc}}

\subsection{Calculation Method\label{sec:calcmethod}}

The purpose of this work is to extend the notion of a mixing parameter~\cite{Johnson1988} to the initialization of the charge density and to bolster the validity of using Janak's theorem~\cite{Janak1978} to calculate the zero-phonon line (ZPL) by demonstrating success for the case of the singly negatively charged calcium vacancy in SiS$_2$, which we investigate as a potential qubit candidate. Our approach builds upon previous work employing Janak's theorem to calculate the electronic properties of excited states of atoms, molecules, and solids.~\cite{PhysRev.184.672,PhysRevA.32.720} In this work, we introduce an innovation by employing Janak's theorem for the calculation of defect properties and additionally providing a lowest order estimate of the error associated with using the theorem for integral change in the occupation of the single-particle states. We argue that the calculations in excess of the ground state calculation in the $\Delta$SCF method are not always necessary, which directly follows from Janak's theorem in the limit where the change in occupation is infinitesimal. To motivate the argument, we compute the error between the $\Delta$SCF method and the approach of using Janak's theorem to lowest order under a change in occupation.

We consider the operator from the single-particle equations in DFT,
\begin{equation}
\mathcal{O}(n(\mathbf{r})) = -\frac{\hbar^2}{2m_e}\nabla^2_{\mathbf{r}}+V(\mathbf{r})+e^2\int\frac{n(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}{\rm d}\mathbf{r}'+\frac{\delta E_{\rm xc}[n(\mathbf{r})]}{\delta n(\mathbf{r}')},
\label{eq:op}
\end{equation}
where $n(\mathbf{r})$ is the particle number density, the first term represents the kinetic energy of noninteracting quasiparticles with electron mass, the second term represents the external potential, the third term represents the Hartree potential for the Coulomb interaction between the quasiparticles and the last term represents the exchange-correlation potential recapturing the fermionic and many-body nature of electronic interactions.

The single-particle equation for the single-particle state $\phi^{(n)}_i$ with eigenvalue $\epsilon^{(n)}_i$ then reads,
\begin{equation}
\mathcal{O}(n(\mathbf{r}))\phi^{(n)}_i = \epsilon^{(n)}_i\phi^{(n)}_i,
\end{equation}
Under a change in occupation, let $n'(\mathbf{r})$ be the new number density and let $\phi^{(n')}_i$ be the new $i^{\rm th}$ single-particle state such that,
\begin{equation}
\mathcal{O}(n'(\mathbf{r}))\phi^{(n')}_i = \epsilon^{(n')}_i\phi^{(n')}_i,
\end{equation}
is satisfied.

The equation for the total energy from DFT is,~\cite{Kohn1965}
\begin{equation}
E = \sum_i^N\epsilon^{(n)}_i -\frac{1}{2}\int\int\frac{n(\mathbf{r})n(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}{\rm d}\mathbf{r}{\rm d}\mathbf{r}'+\int n(\mathbf{r})\left[\epsilon_{\rm xc}[n(\mathbf{r})]-\mu_{\rm xc}[n(\mathbf{r})]\right]d\mathbf{r}.
\end{equation}
If we let,
\begin{equation}
 F[n(\mathbf{r})] = -\frac{1}{2}\int\int\frac{n(\mathbf{r})n(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}{\rm d}\mathbf{r}{\rm d}\mathbf{r}'+\int n(\mathbf{r})\left[\epsilon_{\rm xc}[n(\mathbf{r})]-\mu_{\rm xc}[n(\mathbf{r})]\right]d\mathbf{r},
 \end{equation}
 and
 \begin{equation}
 \Delta \epsilon = \sum_{i=1}^{N+1}(\epsilon^{(n')}_i-\epsilon^{(n)}_i) -(\epsilon^{(n')}_N-\epsilon^{(n)}_N)
 \end{equation}
 then we can approximate the $\Delta$SCF result as,
 \begin{equation}
E_{\rm \Delta SCF} \approx \epsilon^{(n)}_{N+1}-\epsilon^{(n)}_N+ \Delta\epsilon+F[n'(\mathbf{r})]-F[n(\mathbf{r})].
\end{equation}
Therefore, if we only take the difference in ground state eigenvalues, the associated error is,
\begin{equation}
\Delta E_{\rm \Delta SCF} \approx \Delta\epsilon+F[n'(\mathbf{r})]-F[n(\mathbf{r})].
\end{equation}

We can choose to obtain $n'(\mathbf{r})$ at any arbitrary iteration in the full constrained-occupation calculation for the excited state. Therefore, in the limit where $n'(\mathbf{r})$ is obtained at the final iteration of the full constrained-occupation calculation for the excited state, the error becomes exact. In performing the full constrained-occupation calculation for the excited state, we have also explored initializing the charge density by mixing the HOMO and LUMO by varying amounts.

To obtain the defect levels and total energies, we performed first-principles DFT calculations for the various defects in $4H$-SiC and SiS$_2$ using the VASP code~\cite{Kresse1,Kresse2,Kresse3} and the QUANTUM ESPRESSO code~\cite{Giannozzi_2009,Giannozzi_2017} for $\Delta$SCF calculations. In VASP, atomic structures were first converged using the generalized gradient approximation (GGA) for the exchange-correlation energy of electrons, as parametrized by Perdew, Burke and Erzenhof (PBE)~\cite{Perdew2} and then, for the calculation of defect levels of the uncompressed structure, using the screened hybrid functional of Heyd, Scuseria and Ernzerhof (HSE) with the original parameters (0.2 \AA$^{-1}$ for screening and 25\% for mixing).~\cite{Heyd,Krukau} The terms in $F[n(\mathbf{r})]$ were obtained for the ground state and for the state with the changed occupation, where a non-SCF calculation was performed until convergence keeping the charge density fixed in the latter case. The eigenvalues, $\epsilon^{(n')}_i$, were also provided by the non-SCF calculation. Code from work by Feenstra \textit{et al.}~\cite{Feenstra2013} was used to change orbital occupations. In QUANTUM ESPRESSO, we performed $\Delta$SCF calculations to investigate the SiS$_2$ system using PAW pseudopotentials~\cite{Kresse3} with a 108-atom supercell with gamma-point integration. Modified source code was used to alter the charge density in QUANTUM ESPRESSO. The different values of strain investigated in this work were for biaxial strain. For strain values other than the uncompressed structure, the lattice parameters along the strain directions were rescaled from the lattice parameters of the uncompressed structure. The out-of-plane lattice parameter was held fixed. Constant volume relaxations were then performed for the strained structures as described above. Additional computational details may be found in the Appendix.

For first-principles phonon calculations, performed using Phonopy,~\cite{Togo} the atomic positions in the stoichiometric conventional unit cell were relaxed until the magnitude of the Hellmann-Feynman forces was smaller than $10^{-6}$ eV$\cdot$\AA$^{-1}$ with a Monkhorst-Pack grid of $6\times6\times2$ k-points. Supercells containing 108 atoms ($3\times3\times1$ multiple of the conventional unit cell) with appropriately scaled k-point grids and a cutoff energy of 500 eV where then constructed from the stoichiometric conventional unit cell and used to obtain the force constants to compute the phonon dispersion.

\subsection{Formation Energies}

The formation energies of the calcium vacancy in SiS$_2$ in various charge states were calculated according to the formula,~\cite{zhang1991chemical, RevModPhys.86.253}
\begin{equation}
\label{eq:form_eq}
E_f(q) = E_{\text{def}}(q) - E_0 - \sum_i\mu_in_i + q(E_{\text{VBM}} +E_{\text{F}}) + E_{\text{corr}}(q)
\end{equation}
where $q$ denotes the charge state, $E_{\text{def}}(q)$ is the total energy for the defect supercell with charge state $q$, $E_0$ is the total energy for the stoichiometric neutral supercell, $\mu_i$ is the chemical potential of atom $i$, $n_i$ is a positive (negative) integer representing the number of atoms added (removed) from the system relative to the stoichiometric cell, $E_{\text{VBM}}$ is the absolute position of the valence band maximum, $E_{\text{F}}$ is the position of the Fermi level with respect to the valence band maximum (generally treated as a parameter), and $E_{\text{corr}}(q)$ is a correction term to account for the finite size of the supercell when performing calculations for charged defects.~\cite{Vinichenko} This correction term does not simply treat the charged defect as a point charge, but rather considers the extended charge distribution. The chemical potentials of all the reference elements used in our calculations are listed as follows as a function of their crystal structure and total energy per atom: Si (diamond structure, $-5.42$ eV/atom); S (the total energy of a gas phase S$_8$ molecule was calculated and the sublimation enthalpy was then subtracted,~\cite{Kuate2016,Steudel} $-4.20$ eV/atom); and Ca (face-centered cubic structure, $-2.00$ eV/atom). Consideration of Si-rich or S-rich preparation conditions was made, similar to previous work.~\cite{Kuate}

\subsection{Material Screening Approach}
In order to leverage the vast body of existing knowledge in the field, we utilized electronic databases of two-dimensional materials known as the 2D Materials Encyclopedia~\cite{2Dmatpedia} in a top-down filtering approach and the Computational 2D Materials Database (C2DB)~\cite{Haastrup2018,gjerding2021recent} to find the most promising host materials. The entries from the databases were first selected based on the presence of centrosymmetry in their space groups, yielding 4,822 potential candidates with duplicates discarded. Fulfillment of this criterion allows for the possibility of creating in the host a defect without a dipole moment that could interact with local fields. The following filters were then applied (the number of candidates at each filtration step, discarding duplicates, is listed in brackets),
\begin{enumerate}
\item \textit{PBE band gap above 2 eV (985).} The optical transitions of the qubit candidate implanted in the host must not introduce interference from the electronic states of the host.~\cite{Weber8513} The value is taken to be at least as large as the energy of the optical transitions of well-known point defects such as the N$V^-$ center in diamond and singly negatively charged silicon monovacancies in $4H$-SiC. Given that the PBE~\cite{Perdew2} functional generally underestimates the value of the band gap, such a value allows for the defect levels for transitions similar in energy to those of negatively charged silicon monovacancies in $4H$-SiC and the N$V^-$ center in diamond to be separated from the band edges.
\item \textit{Exfoliation energy below 80 meV/atom (441).} The material must be easy to exfoliate and therefore to fabricate. The value is taken to be commensurate with exfoliation energies of common layered materials such as MoS$_2$ with space group P$\bar{6}$m2.
\item \textit{Does not contain an atom with a nuclear spin (6)}. Decoherence caused by the interaction of nuclear spins with the electronic ones of the qubit candidate must be minimized. We set the cutoff at $< 5\%$ natural abundance of isotopes containing a nuclear spin using as reference $4H$-SiC, which is the less restrictive among the two commonly studied systems of diamond and $4H$-SiC with defects showing long coherence times at room temperature.
\end{enumerate}
The filtering process yielded SiS$_2$ with space group P$2_1$/c as the host candidate with the lowest decomposition energy. The other host candidates left at the final stage of filtration were CaO with space group P4/nmm, SO$_2$ with space group P$2_1$/c, CO with space group Cmme, Si$_2$O$_3$ with space group P$2_1$/m, and S$_4$O$_9$ with space group P$\bar{3}$. To our knowledge, none of these compounds have been used as qubit hosts. We note that similar criteria have been used in the work of Ferrenti \textit{et al.},~\cite{Ferrenti2020} though in that work the condition regarding the natural abundance of isotopes with a nuclear spin was relaxed to a cutoff of $<50\%$ and bulk materials were considered, yielding more known candidate hosts. Given the possibility of employing isotope purification, we invite researchers wishing to follow up on our work to relax the cutoff if the above candidates do not prove readily amenable to synthesis.

Working from the literature on defects in diamond and $4H$-SiC, we investigated vacancies, germanium-vacancy complexes and lead-vacancy complexes as possible qubit candidates in the SiS$_2$ host. However, these all break inversion symmetry when they are relaxed with spin polarization. On the other hand, the singly negatively charged calcium-vacancy complex is able to preserve inversion symmetry upon relaxation with spin polarization, due in part to the large size of the calcium atom.

\section{RESULTS AND DISCUSSION \label{sec:res}}

\subsection{Silicon Monovacancies in $4H$-SiC\label{sec:VSi}}
We begin by demonstrating the accuracy of Janak's theorem for silicon monovacancies in $4H$-SiC. This system is an example of a success of the theorem though for another prominent system, the N$V^-$ in diamond, Janak's theorem alone is not sufficient to produce good agreement between experiment and theory. The lack of agreement in that case can be explained by the significant difference between the ground and excited state wavefunctions~\cite{Gali} at a location of significant change in the external or ionic potential (at the N atom in that case), which, from the single-particle equations, will in general lead to larger energy changes under a change in occupation compared to differences elsewhere. We can motivate this argument explicitly using what is referred to as ``effective-mass theory'',~\cite{Kaxiras,Lannoo1992} where the potential introduced by the impurity is considered to be a small perturbation to the crystal potential. Let $V^{\rm defect}(\mathbf{r})$ be the change in the external or ionic potential of the crystal due to the presence of the defect and expand the defect wavefunction of interest, $\phi^{\rm defect}(\mathbf{r})$, in a basis of the eigenfunctions, $\psi_{\mathbf{k}}(\mathbf{r})$, of the perfect crystal,
\begin{equation}
\phi^{\rm defect}(\mathbf{r}) = \sum_{\mathbf{k}}\beta_{\mathbf{k}}\psi_{\mathbf{k}}(\mathbf{r}).
\end{equation}
The wavefunction will obey the single-particle equation,
\begin{equation}
\left[\mathcal{H}^{\rm crystal} + V^{\rm defect}(\mathbf{r})\right] \phi^{\rm defect}(\mathbf{r}) = \epsilon^{\rm defect} \phi^{\rm defect}(\mathbf{r}),
\end{equation}
where $\mathcal{H}^{\rm crystal}$ is the single-particle hamiltonian for the perfect crystal so that,
\begin{equation}
\mathcal{H}^{\rm crystal} \psi_{\mathbf{k}}(\mathbf{r}) = \epsilon_{\mathbf{k}} \psi_{\mathbf{k}}(\mathbf{r}).
\end{equation}
Using the orthogonality of the crystal wavefunctions and the expansion for $\phi^{\rm defect}(\mathbf{r})$ we then have,
\begin{equation}
\epsilon_{\mathbf{k}}\beta_{\mathbf{k}} + \sum_{\mathbf{k}'}\left<\psi_{\mathbf{k}}\right|V^{\rm defect}\left|\psi_{\mathbf{k}'}\right>\beta_{\mathbf{k}'} = \epsilon^{\rm defect}\beta_{\mathbf{k}}.
\end{equation}
Thus, for perturbations to the defect wavefunction that equally change the projection of the defect wavefunction onto some perfect crystal eigenfunction (ie. that result in the same $\beta_{\mathbf{k}}$ for some $\mathbf{k}$), the perturbation that will in general cause the greatest change in energy is the one where the other coefficients are such that the defect wavefunction changes most at the location of the defect potential. As we will see in Section \ref{sec:CaV}, the partial charge densities associated with the orbitals changing occupation can be a quick way to identify when Janak's theorem is likely to fail.

For the lattice parameters of the stoichiometric hexagonal unit cell of $4H$-SiC using the HSE06 functional, we find $a = 3.08$~\AA~and $c = 10.04$~\AA. These values are in good agreement with experimental values~\cite{Levinshtein} of $a = 3.07$~\AA~and $c = 10.05$~\AA~and other theoretical values~\cite{Yan2020} of $a = 3.07$~\AA~and $c = 10.05$~\AA. The structure of the $V_{\rm Si}$ is shown in Fig. \ref{fig:structVSi}. We observe that the expected $C_{3v}$ symmetry is weakly broken from the distances displayed in Table \ref{tab:VSidists}.

\begin{figure}[ht!]
\centering
\includegraphics[width=0.25\textwidth]{structure_VSi.pdf}
\caption{ Shown above is the local structure of a $V_{\rm Si}$ in $4H$-SiC, which can be located at either of the two inequivalent $k$ or $h$ sites. Nearest-neighbor carbon atoms to the vacancy are shown in brown and labeled with numbers and the silicon vacancy is shown as a blue dashed circle. The $V_{\rm Si}$ is displayed along the $[001]$ direction.}
\label{fig:structVSi}
\end{figure}

\vspace{100pt}
\begin{table}[ht!]
\caption{Distances in \AA,~$d(X_1-X_2)$, between carbon atoms numbered in Fig. \ref{fig:structVSi} and the silicon vacancy for $V_{\rm Si}$ at the $k$ and $h$ sites.}
\centering
\vspace{1 mm}
\begin{tabular}{|c|c|c|c|c|}
\hline\hline
site &  $d(V_{\rm Si}-{\rm C}_1)$ & $d(V_{\rm Si}-{\rm C}_2)$ & $d(V_{\rm Si}-{\rm C}_3)$ & $d(V_{\rm Si}-{\rm C}_4)$ \\
\hline
 $k$ &      $2.05$  &   $2.06$ & $2.04$ & $2.03$ \\
 $h$ &      $2.04$  &   $2.04$ & $2.03$ & $2.02$  \\
\hline
\hline
\end{tabular}
\label{tab:VSidists}
\end{table}

From the work of Soykal \textit{et al.},~\cite{Soykal} we know that for a spin-polarized system the approach of considering holes is equivalent to the approach of considering electrons. Defect levels calculated using the HSE06 functional for the singly negatively charged silicon monovacancy $V_{\rm Si}^-$ in $4H$-SiC with $S = 3/2$ at the two inequivalent $h$ and $k$ sites are illustrated schematically in Fig. \ref{fig:eigenVSi} and are provided in Tables \ref{tab:HSEvalsVSik} and \ref{tab:HSEvalsVSih}. We note that the zero of energy is arbitrary as it depends on the pseudopotential, but, as we consider differences in eigenvalues, the value of the zero is not important. From the work of Soykal \textit{et al.},~\cite{Soykal} in the hole picture the ground state manifold is composed of the states, $\left||ue_xe_y+i\bar{u}\bar{e}_x\bar{e}_y\right>/\sqrt{2}$, $\left||ue_xe_y-i\bar{u}\bar{e}_x\bar{e}_y\right>/\sqrt{2}$, $\left||ue_x\bar{e}_y+u\bar{e}_xe_y+\bar{u}e_xe_y\right>/\sqrt{3}$, $\left||\bar{u}\bar{e}_xe_y+\bar{u}e_x\bar{e}_y+u\bar{e}_x\bar{e}_y\right>/\sqrt{3}$, while the excited state manifold in the hole picture is composed of the states, $\left||ve_xe_y+i\bar{v}\bar{e}_x\bar{e}_y\right>/\sqrt{2}$, $\left||ve_xe_y-i\bar{v}\bar{e}_x\bar{e}_y\right>/\sqrt{2}$, $\left||ve_x\bar{e}_y+v\bar{e}_xe_y+\bar{v}e_xe_y\right>/\sqrt{3}$, \\$\left||\bar{v}\bar{e}_xe_y+\bar{v}e_x\bar{e}_y+v\bar{e}_x\bar{e}_y\right>/\sqrt{3}$. Above, $u$ and $v$ are single-particle orbitals transforming as the $A_1$ irreducible representation of the $C_{3v}$ point group, while $e_x$ and $e_y$ transform as the $x$ and $y$ components of the two-dimensional $E$ irreducible representation of the $C_{3v}$ point group. The overbar denotes the minority spin state. To calculate the ZPL we take the lowest energy hole states, which from Tables \ref{tab:HSEvalsVSik} and \ref{tab:HSEvalsVSih} we see must be the $\left||\bar{u}\bar{e}_xe_y+\bar{u}e_x\bar{e}_y+u\bar{e}_x\bar{e}_y\right>/\sqrt{3}$ (excited) and $\left||\bar{v}\bar{e}_xe_y+\bar{v}e_x\bar{e}_y+v\bar{e}_x\bar{e}_y\right>/\sqrt{3}$ (ground) states. Explicitly, in the following we denote the eigenvalue for the majority spin single-particle orbital $\gamma$ as $\epsilon_\gamma$ and for the minority spin single-particle orbital $\bar{\gamma}$ as $\epsilon_{\bar{\gamma}}$. Then, the ZPL we obtain for the $k$ site $V_{\rm Si}^-$ is $(2(\epsilon_{\bar{v}}-\epsilon_{\bar{u}})/3+(\epsilon_{v}-\epsilon_{u})/3) = 1.34$~eV (the remaining terms cancel), while the ZPL we obtain for the $h$ site $V_{\rm Si}^-$ is $(2(\epsilon_{\bar{v}}-\epsilon_{\bar{u}})/3+(\epsilon_{v}-\epsilon_{u})/3) = 1.43$~eV (where again the remaining terms cancel) in excellent agreement with experimental values of 1.35~eV and 1.44~eV, respectively.~\cite{Bracher} The slight underestimation of the ZPL values may therefore be due to the slightly smaller HSE06 band gap (3.18~eV compared to about 3.2~eV for experiment~\cite{Paufler}). We note, however, that other theoretical calculations yield 1.44~eV for the $k$ site and 1.54~eV for the $h$ site using the HSE06 functional.~\cite{Ivady2017} We believe error compensation in taking the difference of many eigenvalues may be causing the greater accuracy of our approximation to the $\Delta$SCF method using Janak's theorem than the $\Delta$SCF method itself, though later theoretical work shows better agreement with our work and with experiment.~\cite{Udvarhelyi2020}

We note that, as the many-body states from the work of Soykal \textit{et al.}~\cite{Soykal} have been defined as Slater determinants of products of single-particle orbitals, the only requirement for the correctness of the calculation outlined above is that these single-particle orbitals be orthogonal, which is consistent with the orbitals constructed in the work of Soykal \textit{et al.}~\cite{Soykal} Though after performing structural relaxation with the HSE06 functional the expected $C_{3v}$ symmetry is weakly broken, we have used the assignment and energetic ordering of the orbitals from the work of Soykal \textit{et al.}~\cite{Soykal} as such an ordering and assignment is consistent with obtaining the lowest energy spin-conserving excitation of the $S = 3/2$ states that our calculations produce. Ultimately, these results demonstrate that in cases where the wavefunctions for the states that will change occupation do not appreciably differ at locations of significant change in the ionic potential and where the change in the ionic potential can be viewed as a small perturbation, there is no need to calculate excited state properties to obtain the ZPL as values derived entirely from the ground state calculation can provide excellent agreement with experiment. These conditions are satisfied for the single vacancy as none of the defect wavefunctions are appreciable at the location of the removal of the atom and a single vacancy is a perturbation of less than $1\%$ of the integrated ionic potential for the supercell size.

\begin{figure}[ht!]
\centering
\includegraphics[width=0.3\textwidth]{VSieigs3.pdf}
\caption{Schematic of the majority (without overbar) and minority (with overbar) spin energy levels calculated in the ground state using the HSE06 functional for the $V_{\rm Si}^-(k)$ and $V_{\rm Si}^-(h)$ point defects in $4H$-SiC. The conduction band is indicated in blue and the valence band in red. Single-particle orbitals transforming as the $A_1$ irreducible representation of the $C_{3v}$ point group are represented by $u$ and $v$, while $e_x$ and $e_y$ transform as the $x$ and $y$ components of the two-dimensional $E$ irreducible representation of the $C_{3v}$ point group. The vertical axis of the figure is not drawn to scale.}
\label{fig:eigenVSi}
\end{figure}

\begin{table}[ht!]
\caption{DFT eigenvalues in eV calculated in the ground state using the HSE06 functional for the hole single-particle states of the $V_{\rm Si}^-(k)$ point defect in $4H$-SiC.}
\centering
\vspace{1 mm}
\begin{tabular}{|c|cc|}
\hline\hline
single-particle state &  majority spin & minority spin \\
\hline
 $u$ &      $7.553$  &   $8.116$  \\
 $v$ &      $7.773$  &  $10.011$ \\
 $e_x$ &      $7.828$  &  $10.152$ \\
 $e_y$ &      $7.855$  &  $10.225$ \\
\hline
\hline
\end{tabular}
\label{tab:HSEvalsVSik}
\end{table}

\begin{table}[ht!]
\caption{ DFT eigenvalues in eV calculated in the ground state using the HSE06 functional for the hole single-particle states of the $V_{\rm Si}^-(h)$ point defect in $4H$-SiC.}
\centering
\vspace{1 mm}
\begin{tabular}{|c|cc|}
\hline\hline
single-particle state &  majority spin & minority spin \\
\hline
 $u$ &    $7.551$  &  $8.119$  \\
 $v$ &    $7.800$  & $10.137$ \\
 $e_x$ &    $7.871$  & $10.227$ \\
 $e_y$ &    $7.906$  & $10.289$ \\
\hline
\hline
\end{tabular}
\label{tab:HSEvalsVSih}
\end{table}

\vspace{110pt}
Based on the results of the calculations for silicon monovacancies in $4H$-SiC outlined above, Janak's theorem shows promise for rapid estimation of ZPL values. Indeed, we will see below that the accuracy of Janak's theorem is maintained for the strain ($\varepsilon$) values with the two lowest total energies, namely the uncompressed and the $\varepsilon = +2\%$ structures, for the singly negatively charged calcium vacancy in SiS$_2$ and that the error calculations we have outlined in Section \ref{sec:calcmethod} provide a clear indication of when Janak's theorem fails.

\subsection{Singly Negatively Charged Calcium Vacancy in SiS$_2$\label{sec:CaV}}
As alluded to above, we now turn to demonstrating the continued accuracy of Janak's theorem for the $\varepsilon$ values with the two lowest total energies for the singly negatively charged calcium vacancy in SiS$_2$. We additionally show that the associated error consistently identifies the larger discrepancies between the results of the theorem and the results of $\Delta$SCF calculations. We also introduce results showing faster convergence of excited state SCF calculations when the charge density is initialized by mixing the HOMO and LUMO compared to when it is initialized from a superposition of atomic charges. For the lattice parameters of the stoichiometric unit cell SiS$_2$ with space group P$2_1$/c using the PBE functional we find $a = 5.93$~\AA, $b = 8.13$~\AA, $\alpha = 90^{\circ}$, $\beta = 102.57^\circ$, $\gamma = 90^\circ$ and a monolayer thickness of $3.65$~\AA~with a vacuum of $16.94$~\AA. Using the HSE06 functional we find $a = 5.88$~\AA, $b = 8.04$~\AA, $\alpha = 90^{\circ}$, $\beta = 103.11^\circ$, $\gamma = 90^\circ$ and a monolayer thickness of $3.59$~\AA~with a vacuum of $16.67$~\AA. The lattice parameters with strain are included in Table \ref{tab:SiS2latparam}.
\begin{table}[ht!]
\caption{PBE lattice constants in~\AA,~and corresponding angles, with strain. The monolayer thickness was not computed as the lattice parameters were directly applied to the defect-containing supercell structures, which were then relaxed at constant volume.}
\centering
\vspace{1 mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline\hline
 $\varepsilon$ &  $a$ & $b$ & $c$ & $\alpha$ & $\beta$ & $\gamma$ \\
\hline
$+5\%$ & $6.22$ & $8.53$ & $21.09$ & $90^\circ$ & $102.57^\circ$ & $90^\circ$ \\
$+2\%$ & $6.05$ & $8.29$ & $21.09$ & $90^\circ$ & $102.57^\circ$ & $90^\circ$  \\
$-2\%$ & $5.81$ & $7.96$ & $21.09$ & $90^\circ$ & $102.57^\circ$ & $90^\circ$ \\
$-5\%$ & $5.63$ & $7.72$ & $21.09$ & $90^\circ$ & $102.57^\circ$ & $90^\circ$  \\
\hline
\hline
\end{tabular}
\label{tab:SiS2latparam}
\end{table}
Structures for the singly negatively charged calcium vacancy for five strain values are found in Fig. \ref{fig:structCaV}. The values in Table \ref{tab:CaVdists} show that compressive strain tends to destroy the $C_i$ symmetry present in the uncompressed, $\varepsilon = +2\%$ and $\varepsilon =+5\%$ structures, due to buckling of the system.

\begin{figure}[ht!]
\centering
\includegraphics[width=0.9\textwidth]{structures.pdf}
\caption{Shown above is the structure of a singly negatively charged calcium vacancy using the PBE functional for $\varepsilon = \pm2, \pm5$ and the uncompressed structure. Silicon atoms are in blue, sulfur atoms are in yellow and the calcium atom is in cyan. Silicon and sulfur vacancies are shown as blue and yellow dashed circles, respectively, for the uncompressed structure. The sulfur atoms closest to the calcium atom in the uncompressed structure are numbered as 1-4.}
\label{fig:structCaV}
\end{figure}

\begin{table}[ht!]
\caption{PBE distances in \AA,~$d(X_1-X_2)$, between sulfur atoms numbered in Fig. \ref{fig:structCaV} and the calcium atom for $\varepsilon = \pm2, \pm5$ and the uncompressed structure as well as the smallest angles, $\theta(X_1-X_2-X_3)$, between the bonds joining the calcium atom to opposite sulfur atoms. HSE06 values are in brackets.}
\centering
\vspace{1 mm}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline\hline
 $\varepsilon$ &  $d({\rm Ca}-{\rm S}_1)$ & $d({\rm Ca}-{\rm S}_2)$ & $d({\rm Ca}-{\rm S}_3)$ & $d({\rm Ca}-{\rm S}_4)$ & $\theta({\rm S}_1-{\rm Ca}-{\rm S}_2)$ & $\theta({\rm S}_3-{\rm Ca}-{\rm S}_4)$ \\
\hline
$+5\%$ & $2.76$ & $2.76$ & $2.83$ & $2.83$ & $179.99^\circ$ & $179.99^\circ$ \\
$+2\%$ & $2.75$ & $2.75$ & $2.80$ & $2.80$ & $179.99^\circ$ & $179.98^\circ$  \\
$0\%$ & $2.76$ ($2.75$) & $2.76$ $(2.75)$ & $2.79$ ($2.79$) & $2.79$ ($2.79$) & $179.98^\circ$ ($179.93^\circ$) & $179.97^\circ$  ($179.89^\circ$)\\
$-2\%$ & $2.76$ & $2.76$ & $2.79$ & $2.79$ & $179.92^\circ$ & $179.88^\circ$  \\
$-5\%$ & $2.97$ & $2.70$ & $3.03$ & $2.72$ & $150.28^\circ$ & $138.39^\circ$  \\

\hline
\hline
\end{tabular}
\label{tab:CaVdists}
\end{table}

The approximation to the ZPL using Janak's theorem and using the $\Delta$SCF method for different in-plane $\varepsilon$ can be found in Fig. \ref{fig:ZPLCavals}, where the corresponding eigenvalues can be found in Table \ref{tab:HSEPBEvalsSiS2}. The difference between the LUMO and the HOMO eigenvalues is used in the application of Janak's theorem. We observe little variation of the approximation to the ZPL with tensile $\varepsilon$, but more variation with compressive $\varepsilon$ which can also distort the structure. Since, due to Poisson's ratio, the application of pressure should lead to tensile in-plane $\varepsilon$, we do not expect the ZPL value to change significantly from what we have predicted in experimentally realizable structures. The error calculations perform best for the uncompressed (error of $0.0128$~eV) and $\varepsilon = +2\%$ (error of $0.0249$~eV) structures, which have the smallest energy differences, obtaining the correct sign of the error as well, but fail to capture the correct sign and are much too large for the remaining $\varepsilon$ values (errors of $-17.4096$~eV, $-0.3898$~eV and $-0.6123$~eV for $\varepsilon = -5\%, -2\%$ and $+5\%$, respectively). We note that for $\varepsilon = -5\%$, the Ca atom departed from the position that preserved inversion symmetry and that strain value was therefore not similar to the other structures. The larger errors nonetheless correctly indicate more significant disagreement between the $\Delta$SCF result and the result from using Janak's theorem. We note that based on a statistical learning based prediction of the bulk modulus,~\cite{deJong2016,Evers2015} the structures that should have the two lowest total energies given the requirement of an applied pressure of 2.8 GPa to 3.5 GPa are the uncompressed and $\varepsilon = +2\%$ structures, which is verified by QUANTUM ESPRESSO total energy calculations. Naturally, the uncompressed structure has the lower energy of the two structures.

As shown in Fig. \ref{fig:beta_conv}, these structures with the two lowest total energies also responded best to performing the excited state calculation by replacing the initialization of the charge density from a superposition of atomic charges with the charge density of the ground state calculation, where the HOMO contribution is modified according to,
\begin{equation}
\label{eq:betaHL}
\left|{\rm HOMO}\right> \rightarrow \beta\left|{\rm HOMO}\right>+(1-\beta)\left|{\rm LUMO}\right>,
\end{equation}
$0\leq \beta \leq 1$. We note that a ground state relaxation is needed for both the standard method or default initialization and the approach of initializing the charge density from a mixture of the HOMO and LUMO orbitals for the excited state SCF calculation. This prerequisite calculation therefore does not change the net increase or decrease in the number of iterations for one method over another. Given that there may be instances where a ground state relaxation is not performed, we report the number of iterations required to converge the ground state SCF calculation for completeness. For $\varepsilon = -5\%, -2\%, +2\%, +5\%$ and the uncompressed structure, the number of iterations required were 29, 35, 34, 31, and 35, respectively.
The value of the band-gap for the uncompressed structure using the PBE functional is 3.55~eV and the defect has total spin $S = 1/2$. The $\Delta$SCF calculations and the differences in ground state eigenvalues were from QUANTUM ESPRESSO.

\begin{figure}[ht!]
\centering
\includegraphics[width=0.55\textwidth]{ZPL_bars.pdf}
\caption{ZPL values in eV, calculated using Janak's theorem (red) and using the $\Delta$SCF method (green), for the singly negatively charged calcium vacancy in SiS$_2$ for in-plane $\varepsilon = \pm2\%, \pm5\%$ and for the uncompressed structure. The error associated with the result from Janak's theorem was calculated as described in Section \ref{sec:disc}. The $\varepsilon = -5\%$ case caused the Ca atom to depart from the position that preserved inversion symmetry and was therefore not similar to the other structures.}
\label{fig:ZPLCavals}
\end{figure}

\begin{table}[ht!]
\caption{DFT eigenvalues in eV calculated in the ground state of the negatively charged calcium vacancy in SiS$_2$ using the PBE and HSE06 functionals for various $\varepsilon$ values.}
\centering
\vspace{1 mm}
\begin{tabular}{|c|ccccc|c|}
\hline\hline
& \multicolumn{5}{c|}{ PBE}  & HSE06 \\
\hline
  $\varepsilon$ & $-5$\% & $-2$\% & uncompressed & $+2$\% & $+5$\% & uncompressed \\
\hline
 HOMO &  $-2.2137$  & $-2.8524$  & $-3.1014$ & $-3.2690$ & $-3.4880$ & $-3.1995$   \\
 LUMO  &  $-1.7557$  & $-2.7062$  & $-2.9211$ & $-3.0855$ & $-3.3124$ & $-2.9319$   \\
\hline
\hline
\end{tabular}
\label{tab:HSEPBEvalsSiS2}
\end{table}

\begin{figure}[ht!]
\centering
\includegraphics[width=0.5\textwidth]{beta_conv2.pdf}
\caption{Change in the number of iterations to convergence using the charge density initialization presented in Eq. (\ref{eq:betaHL}) and using the default initialization in QUANTUM ESPRESSO. The zero represents the default initialization or standard method so that a data point with a value along the vertical axis of $-7$ indicates that the particular run converged in $7$ fewer iterations than the standard method for the appropriate strain value. In-plane $\varepsilon = \pm2\%, \pm5\%$, where $\varepsilon$ denotes strain, and the uncompressed structure were investigated. The maximum number of iterations was set at 100, which was attained for both $\varepsilon = +2\%$ and $\varepsilon = +5\%$ and explains the plateaus. We took $0 \leq \beta \leq 1$ and used increments of 0.01. Data points are indicated by `$\times$' and the dashed lines are a guide for the eye. The number of iterations required to converge the ground state SCF calculation is included in the text for all $\varepsilon$ values for completeness.}
\label{fig:beta_conv}
\end{figure}

Based on the computational efficiency of our method for determining the error associated with using Janak's theorem demonstrated in Fig. \ref{fig:times}, we argue that the approach of using Janak's theorem should be acceptable for screening large numbers of potential point defect single-photon emitter candidates for desired ZPL values. Our approach to the initialization of the charge density also shows promise. A key point is that the search would focus on point defect candidates with total spin $S = 1/2$, otherwise more involved group theoretic arguments would be required to determine the correct many-body hole wavefunction as in the work of Soykal \textit{et al.}~\cite{Soykal} Other specific systems, such as triplet systems,~\cite{smart2020intersystem} will also be the focus of future work. The band-gap for the uncompressed structure using the HSE06 functional is 4.75~eV and we again find that the defect has total spin $S = 1/2$. Using QUANTUM ESPRESSO, the HSE06 calculation yields a $\Delta$SCF result of 0.1705~eV and a difference in ground state eigenvalues of 0.2676~eV. The error calculated as described in Section \ref{sec:disc} yields 0.3008~eV. For hosts where the accuracy of hybrid functionals has been demonstrated, clearly their continued use would be recommended. However, one cannot generally state that hybrid functionals will be more accurate than semilocal ones. Indeed, Johari and Shenoy~\cite{Johari2012} find that their PBE band gap for monolayer MoS$_2$ underestimates the experimental band gap of 1.8~eV by just 0.12~eV, while Kuc \textit{et al.}~\cite{Kuc2011} and Ataca and Ciraci~\cite{Ataca2011} found that the PBE0 and HSE06 hybrid functionals overestimate this band gap by approximately 1 and 0.45~eV, respectively. In our work, we have therefore included both the PBE and HSE06 ZPL transitions for the uncompressed structure to reduce the possibility of missing the transition, noting that the transition energy does not vary significantly with strain for low total energy strain configurations. We recommend that any application of our methods using a semilocal functional also be performed with a hybrid functional.

\begin{figure}[ht!]
\centering
\includegraphics[width=0.55\textwidth]{Time_bars.pdf}
\caption{Runtimes for convergence of the constrained-occupation calculation for the excited state (yellow), for completion of 20 SCF iterations (turquoise) and for calculation of the error associated wth using Janak's theorem using our method (blue). As our error calculations were done using the VASP code and the constrained-occupation calculations for the excited state were done using the QUANTUM ESPRESSO code, we took the time for a single SCF iteration for a given system from VASP and multiplied by the number of iterations required for the constrained-occupation calculation for the excited state in QUANTUM ESPRESSO to converge.}
\label{fig:times}
\end{figure}

The DFT partial charge densities obtained for ground hole (LUMO) and excited hole (HOMO) states and differences (LUMO $-$ HOMO) are shown in Fig. \ref{fig:eigenCa}. The differences show the closeness of certain LUMO-HOMO pairs. More importantly, Fig. \ref{fig:eigenCa} shows that Janak's theorem succeeds when the differences are not appreciable at the location of the greatest change in the ionic potential of the perfect lattice (at the location of the Ca atom) and does worse otherwise. The structure of the defect is four missing atoms, two silicon and two sulfur, replaced by a single calcium atom such that the resulting point defect is inversion-symmetric, as confirmed by the partial charge densities in Fig. \ref{fig:eigenCa}. This inversion symmetry is broken, however, for in-plane $\varepsilon = -5\%$. The lowest formation energy as a function of Fermi level is displayed in Fig. \ref{fig:formCa} for the uncompressed structure. The plot demonstrates that in sulfur-rich preparation conditions the introduction of a calcium vacancy actually stabilizes the SiS$_2$ structure with space group P$2_1$/c. Given the fact that the singly negatively charged calcium vacancy only exists for a very limited range of Fermi level values, it would be necessary to pin the Fermi level at the appropriate value by doping the system. Using a $^{43}$Ca atom instead of a $^{40}$Ca atom could also lead to coupling between the nuclear spin and the electronic spin to implement a long-lived quantum memory realized with the nuclear spin.~\cite{Pfender2017}

\begin{figure}[ht!]
\centering
\includegraphics[width=0.5\textwidth]{orbitals_and_tot_lg.pdf}
\caption{DFT calculated partial charge densities for singly negatively charged calcium vacancy in SiS$_2$ using the PBE functional for in-plane $\varepsilon = \pm2\%, \pm5\%$ and for the uncompressed structure and using the HSE06 functional for the uncompressed structure are displayed corresponding to: the LUMO (left), the HOMO (middle) and the difference between the total charge densities with and without exchanging the occupations of the HOMO and LUMO (right). The value of $\varepsilon$ is indicated to the left of the respective panels. Sulfur atoms are in yellow, silicon atoms are in blue and the calcium atom is in cyan. Charge accumulation (depletion) is indicated by translucent red (green). The isosurface of charge density is 0.0005 $e/$\AA$^3$ for all plots.}
\label{fig:eigenCa}
\end{figure}

\begin{figure}[ht!]
\centering
\includegraphics[width=0.5\textwidth]{CaV_form.pdf}
\caption{Formation energy as a function of Fermi level for the calcium vacancy for a PBE DFT calculated gap $E_g = 3.55$ eV. The lower line (yellow) corresponds to sulfur-rich preparation conditions while the upper line (blue) corresponds to silicon-rich preparation conditions. The integers between the lines indicate the most stable charge state of the defect at the corresponding values of the Fermi level, with each region bounded by vertical dotted lines. Calculations were done for the uncompressed structure.}
\label{fig:formCa}
\end{figure}

We now turn to addressing the possibility of non-radiative decay with phonon emission upon excitation of the singly negatively charged calcium vacancy. Unlike in diamond where optical phonon modes with energies of about 0.17~eV can be found,~\cite{Sato2020} the energies of the phonon modes in SiS$_2$ with space group P$2_1$/c are all below 20 THz as shown in Fig. \ref{fig:phondisp}, corresponding to an energy of about 0.08 eV. This result indicates that non-radiative decay would require multi-phonon processes, which would be rare.

 We note that unlike diamond which is routinely synthesized at high pressure but shows no dynamical instability at ambient condition, our calculations show (see Fig. \ref{fig:phondisp}) that SiS$_2$ with space group P$2_1$/c does exhibit a dynamical instability at ambient condition and would therefore need to be maintained at the 2.8 GPa required for its synthesis to be stable. We would like to make the following points: Experiments at high pressure ($> 10^2$~GPa) are routinely performed. In particular, the discovery of room-temperature superconductors at high pressure (e.g. Snider \textit{et al.}~\cite{Snider2020} at 267 GPa) will undoubtedly propel further technical advances in this area. Experiments of emitters at high pressure (e.g. NV centers in diamond) have also been performed (e.g. Lesik \textit{et al.},~\cite{Lesik2019} Yip \textit{et al.},~\cite{Yip2019} and Hsieh \textit{et al.}~\cite{Hsieh2019}), showing such experiments to be within current feasibility. While it is not ideal that the proposed candidate material requires high pressure, it does not remove one of the key elements of the manuscript $-$ the tools for the rapid screening of defect qubits.

\begin{figure}[ht!]
\centering
\includegraphics[width=0.5\textwidth]{pbands_331_662k.pdf}
\caption{Phonon dispersion for the uncompressed P$2_1$/c structure of SiS$_2$ along the indicated high-symmetry path, as calculated with the phonopy code.~\cite{Togo} The structure does show imaginary frequencies near the $\Gamma$ point and would therefore need to be maintained at the 2.8 GPa required for its synthesis to be stable.~\cite{Plasienka2016} The vectors $\mathbf{a^*}$, $\mathbf{b^*}$ and $\mathbf{c^*}$ are reciprocal to the vectors corresponding to the lattice constants $a$, $b$ and the out-of-plane lattice constant, respectively.}
\label{fig:phondisp}
\end{figure}
In summary, the case of the singly negatively charged calcium vacancy in SiS$_2$ confirms that using Janak's theorem and the associated error estimation for ZPL calculations is quick and clearly signals when the theorem is applicable. We have also shown that the constrained occupation calculation for the excited state can be completed in some cases with fewer SCF iterations by carefully initializing the charge density.

\section{CONCLUSION \label{sec:conc}}
In conclusion, we have developed a novel computational framework for high-throughput screening of the ZPL of atom-like qubits with applications for materials discovery, quantum sensing and quantum computing. We propose the use of Janak's theorem for ZPL calculations and provide a quick new method for estimating the associated error. We also outline a new method for reducing the number of iterations required to converge excited state SCF calculations, compared to the default procedure for carrying out such calculations. Our ZPL results are consistent with previous experimental work where available and suggest that considering only ground state eigenvalues in the hole approach is a computationally efficient way of calculating optical excitation energies of color centers for screening purposes. Our novel method for initializing the charge density to reduce the number of iterations required for excited state calculations also works well for systems with small strain. Finally, we propose the new singly negatively charged calcium vacancy in SiS$_2$, which has the advantage of being inversion-symmetric and of containing a low density of nuclear spins that would lead to poor coherence properties of spins hosted in the material. However, this new proposed system would need to be maintained at the 2.8 GPa required for its synthesis to ensure stability. We also note that our rapid screening methods can create training data for machine learning and serve in technologies for predicting defect structures with desirable properties.

\section*{ ACKNOWLEDGMENTS:}

We thank Efthimios Kaxiras of Harvard University and Georgios A. Tritsaris for useful discussions. R.K.D. gratefully acknowledges financial support from the Princeton Presidential Postdoctoral Research Fellowship. We also acknowledge support by the STC Center for Integrated Quantum Materials, NSF Grant No. DMR-1231319. This work used computational resources of the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant Number ACI-1548562,~\cite{Towns} on Stampede2 at TACC through allocation TG-DMR120073, and of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

\section*{DATA AVAILABILITY:}

The data from DFT calculations that support the findings of this study and the modified source code for our charge density implementation are available from Rodrick Kuate Defo (email: rkuatedefo@princeton.edu) upon reasonable request.

\section*{APPENDIX: COMPUTATIONAL DETAILS}

We provide in this Appendix details of the computations. The atomic positions were relaxed until the magnitude of the Hellmann-Feynman forces was smaller than $10^{-4}$~eV$\cdot$\AA$^{-1}$ on each atom without spin-polarization and subsequently until the magnitude of the Hellmann-Feynman forces was smaller than $10^{-2}$~eV$\cdot$\AA$^{-1}$ on each atom with spin polarization to obtain defect levels and, for the stoichiometric conventional unit cell, the lattice parameters were concurrently relaxed. The wavefunctions were expanded in a plane wave basis with a cutoff energy of 500~eV for all systems and a Monkhorst-Pack grid of $6\times6\times2$ k-points was used for integrations in reciprocal space for SiS$_2$ and a Gamma centered grid of $4\times4\times2$ k-points was used for integrations in reciprocal space of $4H$-SiC. The relaxed lattice parameters of the stoichiometric unit cell were then used for the supercell structures. Formation energies and defect levels were calculated using a supercell of 108 atoms for SiS$_2$ ($3\times3\times1$ multiple of the stoichiometric unit cell) with appropriately scaled k-point grids. For $4H$-SiC, a supercell of 576 atoms ($6\times6\times2$ multiple of the stoichiometric unit cell) was used.