\documentclass[a4paper,oneside]{saip}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% defaults, keep as is, don't add anything that changes the document's defaults.
\usepackage{helvet}          % to have better matching fonts between Latex and Word (Title/Headings)
\usepackage{times}           % to have better matching fonts between Latex and Word (bulk text)
\usepackage{courierten}   % to have better matching fonts between Latex and Word (monospaced)
\usepackage[T1]{fontenc} % needed to deal with these fonts

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% common packages needed, add or comment as needed
\usepackage{graphicx}  %to include figures
\usepackage{float} % for [H] float placement
\usepackage{url}
\usepackage{geometry} %	need: 
\usepackage{amsmath} % for \text in math mode
\geometry{%showframe, 
top=40.0mm, bottom=27.0mm, left=25.0mm, right=25.0mm}
% ------------- pre‑amble (before \begin{document}) -----------------
\makeatletter
\def\IEEEbibitemsep{0pt}%            % keep IEEEtran’s stretch at 0

% kill the list’s own \itemsep / \parskip inside the bibliography
\let\oldthebibliography\thebibliography
\let\endoldthebibliography\endthebibliography
\renewenvironment{thebibliography}[1]{%
  \begin{oldthebibliography}{#1}%
    \setlength{\itemsep}{0pt}%       % ← no space between \bibitems
    \setlength{\parskip}{0pt}%       % ← no extra paragraph gap
}{%
  \end{oldthebibliography}%
}
\makeatother

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% optional packages, you can add others as well for your own needs
%\usepackage{xcolor}
%\usepackage{algorithm2e}
%\usepackage[version=4]{mhchem} %for isotopes and all with \ce{^{15}_{8}O_{7}}
%\usepackage{amsmath}
%\usepackage{amssymb}
%\usepackage{stmaryrd}
%\usepackage{mathrsfs}
%\usepackage{mathabx}
%\usepackage{upgreek}
%\usepackage{comment}
%\usepackage{verbatim}
%\usepackage{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DO NOT USE the hyperref package, your PDF will be included in a bigger document and the links will be broken.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

\title{Assessing the Freysoldt, Neugebauer \& van de Walle (FNV) and Kumagai–Oba (KO) finite-size corrections for Ce-vacancy complexes in diamond}

\author{W. Dednam$^{1}$ and E. B. Lombardi$^{1}$}

\affil{$^1$ Research and Academic Computing: High-Performance Computing (HPC), College of Science, Engineering and Technology, University of South Africa, Tshwane, 0001, South Africa}

%only put the email of the main corresponding author, not of all authors.
\email{\url{dednaw@unisa.ac.za}}

\begin{abstract}
Lanthanide-doped diamond couples the ultra wide band gap, high thermal conductivity and radiation hardness of the host with the rich electron physics of the dopant, promising solid-state qubits, single-photon emitters and spin-memory elements. Supercell density functional theory treatments of the charged Ce-vacancy complexes that underlie these functionalities suffer from spurious image–image interactions and an ill defined electrostatic zero energy; post processing finite-size corrections are therefore mandatory for quantitative defect thermodynamics.

We benchmark the two principal correction schemes -- the reciprocal space potential alignment\\monopole method of Freysoldt, Neugebauer, van de Walle (FNV) and the real space multipole expansion of Kumagai Oba (KO) -- for Ce$_{\rm V2}$, Ce$_{\rm V3}$ and Ce$_{\rm V4}$ in a $216$-atom diamond supercell. Uncorrected neutral formation energies agree with literature to within $0.5$eV for all three and confirm Ce$_{\rm V3}$ as the most stable neutral complex. Introducing positive charge exposes limitations of FNV: once the anisotropic Ce $4f$ charge density departs from the isotropic monopole assumed in that formalism, FNV corrections fail to converge. In contrast, KO, which accounts for higher multipoles, remains numerically stable and delivers consistent corrections.

Our results show that KO is indispensable for heavy-atom defects with non-spherical charge distributions, while FNV is reliable only for nearly isotropic cases. This enables accurate assessment of rare-earth dopants in diamond and other wide-gap semiconductors.
\end{abstract}

\section{Introduction}

% This document has instructions for writing your proceedings using LaTeX.  {\bf Your goal is to create a document in the style of} ``\texttt{001-yourname.pdf}''.  You will start with ``\texttt{001-yourname.tex}'', change the ``001'' in the filename with your contribution number (3 digits with leading zeros if needed) and replace ``yourname'' with either your name or surname: {\bf no special characters, no spaces, no symbols}, only A-Z characters without accents.  \\
Doped diamond complexes are promising for spintronic \cite{Lombardi2003,Zhao2011,Benecha2018,Gulka2021} and quantum optics applications \cite{Magyar2014,Sedov2017,Rembold2020,Sedov2022}. Diamond is an attractive host for these defects because it is a radiation-hard, ultra wide band gap semiconductor capable of hosting deep defect states. For example, nitrogen-doped diamond already exhibited promise towards spintronic applications in theoretical density-funtional theory (DFT) calculations in Ref.~\cite{Lombardi2003}, and the negative nitrogen vacancy cluster has been shown to exhibit spin conherence times at room temperature up to the order of milliseconds and spin-relaxation times beyond that, which is very promising for quantum memory applications \cite{Zhao2011}. Subsequent works have continued to explore such advantageous aspects of this defect center (see Ref.~\cite{Gulka2021} and references therein). 

Lanthanide-doped diamond has also recently been explored for its potential, particularly in quantum optics applications. Several studies have explored Cerium \cite{Tan2020Ce,Tan2021}, Erbium \cite{Cajzl2018,Propst2023}, Europium \cite{Magyar2014,Vanpoucke2019,Tan2020Eu}, Praseodymium \cite{Tan2018,Morris2024} and Ytterbium \cite{Cajzl2018,Babar2024} as dopants in diamond complexes of various forms.  

Density-functional theory (DFT) has served as the workhorse for modelling a broad range of diamond complexes \cite{Lombardi2003,Benecha2018,Tan2018,Vanpoucke2019,Tan2020Ce,Tan2020Eu,Cajzl2022,Morris2024,Ghafari2025,Wang2025}. However, the periodic-boundary-condition (PBC) supercell model used in DFT creates spurious electrostatic interactions between a charged defect and its periodic images, so explicit correction terms are required \cite{Makov1995,Lany2008,Freysoldt2009,Komsa2012,Kumagai2014}. Conventional potential-alignment schemes give only a partial correction \cite{Lany2008}; more sophisticated approaches (e.g. image-charge correction schemes) are needed for reliable formation energies \cite{Freysoldt2009,Komsa2012,Kumagai2014}. In this work we assess two of these correction schemes, the \textit{scalar}  Freysoldt, Neugebauer, van de Walle (FNV) and the \textit{tensor} Kumagai Oba (KO), and apply them to the Ce$_{\rm V2}$, Ce$_{\rm V3}$ and Ce$_{\rm V4}$ complexes of diamond studied in Ref.~\cite{Tan2020Ce}.

\section{Methods}

%\subsection{Formation energy of a charge supercell in DFT} 
% Once renamed, edit the ``\texttt{001-yourname.tex}'' file directly and replace what you see in there with your own information, figures, tables, results, etc. Cite your sources properly using bibtex. For example, you might want to cite this paper \cite{1996ApJ...462..563N} or that paper \cite{2003ApJS..148..175S}.  In certain cases, maybe citing them together would make more sense, like this \cite{1996ApJ...462..563N, 2003ApJS..148..175S}.   All the relevant information for these citations are in the ``\texttt{biblio.bib}'' file already.  This template has all the code / packages needed to generate the bibliography correctly already, it's a working example.  Edit as needed by adding more entries for your sources / citations.

The FNV and KO finite-size corrections outlined below follow their implementation in the \texttt{SPINNEY} post-processing code, with conventions as defined by Ref.~\cite{Arrigoni2021}. The formation energy of a charged defect supercell in DFT is conventionally defined as \cite{Arrigoni2021}:

\begin{equation}
E_{f}\!\bigl[D^{q}\bigr] = 
\left\{ E\!\bigl[D^{q}\bigr] + E_{\rm{corr}}\!\bigl[D^{q}\bigr] \right\}
- E_{P}
-\sum_{i} n_{i}\mu_{i}
+ q\!\left( \varepsilon_{\mathrm{VBM}} + \Delta\varepsilon_{F} \right).
\label{eq:FE}   
\end{equation}

where $E_{f}\!\bigl[D^{q}\bigr]$ is the ground-state DFT total energy of the defect supercell in charge state $q$ and $E_{P}$ is the neutral pristine diamond supercell.
$n_{i}$ is the number of removed ( $< 0$) or added ( $> 0$) atoms  and $\mu_{i}$ is the chemical potential of atom type $i$.
$E_{\rm VBM}$  is the valence-band maximum (VBM) of the pristine supercell. $\Delta\varepsilon_{F}$ is the Fermi level relative to the VBM. %, used to plot charge-transition levels (CTLs) of the defects up to the conduction-band minimum (CBM) of the pristine supercell.
Finally, $E_{\rm{corr}}\!\bigl[D^{q}\bigr]$ consists of two terms: the image-charge correction and potential alignment. 

%\subsection{Why corrections to the formation energy are needed to begin with}

Under PBC, a charged defect interacts with its own periodic images and the neutralising background ($\frac{-q}{\Omega}$, with $\Omega$ the supercell volume). Image-charge corrections remove this artificial long-range interaction, so that $E_{f}\!\bigl[D^{q}\bigr]$ converges toward the dilute-defect limit in the ideal case \cite{Kumagai2014}. A potential alignment correction is also needed since the DFT Poisson solver sets the cell-average potential $V_{{\rm avg},\Omega}$ to zero arbitrarily \cite{Kumagai2014}. The charged cells inherit this potential shift, which, in turn, alters Kohn-Sham eigenvalues and total energies.

% Sometimes what you want to cite isn't an article.  For example, this is how to cite an online instruction manual \cite{echelle}; check the bib file and the references at the end.   Also, this proceedings template is based on the one used by the ``Journal of Physics: Conference Series'' \cite{jpconf}.  It has been modified and expanded for the specific needs of the SAIP Conference Proceedings, all contained in the ``\texttt{saip.cls}'' file.  


%\subsection{The \textit{scalar} FNV correction}

The FNV method defines the correction terms as \cite{Arrigoni2021}: 

\begin{equation}
E_{{\rm corr}}\!\bigl[D^{q}\bigr]
   = -\,E_{\rm model}
     - q\,\bigl.\Delta V_{{\rm model},\,q/b}\bigr|_{\rm{far}}
\label{eq:FNV}   
\end{equation}

where $E_{\rm model}$ is the electrostatic (Madelung) energy of a model charge density, which can be a point charge (PC) or a charge density fitted to a suitable (radially isotropic) function.
 $\Delta V_{{\rm model},\,q/b}\bigr|_{\rm{far}}$ is the value of $\Delta V_{{\rm model},\,q/b}$, the difference between the defect-induced potential $V_{q/b}$ and the and the charge-density model potential $V_{{\rm model,}q}$ far away from the defect inside $\Omega$. For convenience, the overall correction term in Eq.~(\ref{eq:FE}) has traditionally been analyzed in Ewald summation form divided into real and reciprocal space parts \cite{Makov1995,Kumagai2014}. The electrostatic component of this correction can thus be defined as \cite{Kumagai2014,Arrigoni2021}:

\begin{equation}
E_{\rm model}
   = \frac{1}{2\,\varepsilon\,\Omega}
     \sum_{\mathbf G \neq \mathbf 0}
     \frac{4\pi}{|\mathbf G|^{2}}\,
     \bigl|\tilde\rho_{\rm model,d}(\mathbf G)\bigr|^{2}
     \;-\;
     \frac{1}{\pi\,\varepsilon}
     \int_{0}^{\infty}
       \kappa\,
       \bigl|\tilde\rho_{\rm model,d}(\kappa)\bigr|^{2}\,d\kappa ,
\label{eq:FNVmodel}   
\end{equation}

where $ \tilde\rho_{\rm model,d}(\mathbf G)
   = \int_{\Omega}
       \rho_{\rm model,d}(\mathbf r)\,
       e^{-i\mathbf G \cdot \mathbf r}\,d^{3}r$, $\kappa = |\mathbf G|$, and $\varepsilon$ is the dielectric constant of the relaxed defect supercell calculated within the linear response regime in our DFT calculations. The first term of Eq.~(\ref{eq:FNVmodel}) is just the Fourier transform of the electrostatic energy of $\rho_{\rm model,d}$ and the second term corrects for the self-energy of an isolated charge density \cite{Kumagai2014}. The model charge density is defined as \cite{Freysoldt2009,Arrigoni2021}:

\begin{equation}
\rho_{\rm model,d}(\mathbf r)
  = qx\,\frac{e^{-|\mathbf r - \mathbf r_{0}|/\gamma}}{N_{\gamma}}
    + q(1-x)\,\frac{e^{-|\mathbf r - \mathbf r_{0}|^{2}/\beta^{2}}}{N_{\beta}},
\label{eq:FNVrho}   
\end{equation}    

where $N_{\gamma} = 8\pi\,\gamma^{3}$, $N_{\beta}  = \pi^{3/2}\,\beta^{3}$, and $x$, $\beta$ and $\gamma$ are fitting parameters with ${\mathbf r}_0$ is centered on the defect \cite{Freysoldt2009,Komsa2012,Kumagai2014}. Equation~(\ref{eq:FNVrho}) represents a weighted sum of an exponentially decaying radial charge distribution and a localised Gaussian charge distribution. If the weighting factor $x=1$, then the model charge density is completely delocalised, with the opposite case occurring when $x=0$: the charge density is fully localised.  

Importantly, the potential alignment $\Delta V_{{\rm model},\,q/b}\bigr|_{\rm{far}}$ in the FNV correction varies strongly around the defect as the extent of relaxation of atoms around it increases \cite{Komsa2012,Kumagai2014}.
A planar average of this potential is used in the FNV correction \cite{Freysoldt2009,Arrigoni2021}. Hence, relaxation of atoms at the defect gives rise to an anisotropic charge distribution which the scalar FNV correction may not be able to accurately capture \cite{Komsa2012,Kumagai2014}.

% We are using ``IEEE-style'' for our bibliography, from the Institute of Electrical and Electronics Engineers (IEEE) transactions journals and conferences (``\texttt{IEEEtran.bst}'')\footnote{Footnote with link: \url{https://www.ctan.org/tex-archive/macros/latex/contrib/IEEEtran}}.  Sometimes, it might be more appropriate to put a footnote with an link rather than a full reference at the end, as shown here.  \\

%\subsection{The \textit{tensor} KO extension}

To remedy this inherit deficiency, the \textit{tensor} KO extension of the \textit{scalar} FNV correction was introduced \cite{Kumagai2014}:

\begin{equation}
E_{\rm corr}\!\bigl[D^{q}\bigr]
   = -\,E_{\rm PC}
     - q\,\Delta V_{\rm PC,q/b}
\label{eq:KO}   
\end{equation}

\begin{equation}
E_{\rm PC}
   = \frac{1}{2\Omega}\,
     \sum_{\mathbf G\neq\mathbf 0}
     \frac{4\pi}{\mathbf G^{\mathrm T}\!\mathbf{\varepsilon}\,\mathbf G}\,
     \bigl|\tilde\rho_{\rm model,d}(\mathbf G)\bigr|^{2}
   = \frac{2\pi q^{2}}{\Omega}\,
     \sum_{\mathbf G\neq\mathbf 0}
     \frac{1}{\mathbf G^{\mathrm T}\!\mathbf{\varepsilon}\,\mathbf G}
\label{eq:KOmodel}   
\end{equation}

\begin{equation}
\tilde\rho_{\rm model,d}(\mathbf G)
   = \int_{\Omega}
       \rho_{\rm model,d}(\mathbf r)\,
       e^{-i\mathbf G\!\cdot\!\mathbf r}\,d^{3}r
   = \int_{\Omega}
       q\,\delta(\mathbf r-\mathbf r_{0})\,
       e^{-i\mathbf G\!\cdot\!\mathbf r}\,d^{3}r
   = q,
\label{eq:KOrho}   
\end{equation}     

where now $\mathbf \varepsilon$ is the dielectric \textit{tensor} also obtained from a DFT linear response calculation on the relaxed defect supercell. In practice, $\mathbf{\varepsilon}$ is diagonalized and both $\mathbf{\varepsilon}$ and $\mathbf{G}$ are rotated into the principal-axis frame (where $\mathbf{\varepsilon}$ is diagonal), so that $\mathbf{G}^{\mathrm{T}} \mathbf{\varepsilon}\, \mathbf{G} = \sum_i \varepsilon_i\, G_i^2$ \cite{Kumagai2014}. Note that here the constant (and divergent) electrostatic self-interaction term cancels exactly in Eq.~(\ref{eq:KO}) because it appears with the same magnitude but opposite sign in the electrostatic energy $E_{\rm PC}$ and the reference potential $V_{\rm PC}$ used to define the alignment term $\Delta V_{\rm PC,q/b}^{\rm far}$ in the KO scheme \cite{Kumagai2014}. Hence, no explicit self-energy contribution is needed when Eq.~(\ref{eq:KO}) is used in Eq.~(\ref{eq:FE}).

% \label{sec:fig-tab}
% Include your figures and tables and refer to them in the order they will appear, figure 1 before figure 2, etc.  This is how to refer to Table~\ref{tab:mytable} and Figure~\ref{fig-sidebyside}, Figure~\ref{fig-stackup}, and Figure~\ref{fig-onlyone}.  No figure or table should be displayed in your paper without being referred to in the text. \\

%\subsection{Ce-vacancy complexes in diamond}

Figure~\ref{fig1} shows the Ce-vacancy complexes considered in this work: Ce inserted into a divacancy, a trivacancy or at the site of 4 C atoms removed. They have been created from a relaxed $216$-atom diamond supercell with lattice constant $3.574\,\text{\AA}$, in close analogy with the structures employed in Ref.~\cite{Tan2020Ce}. \footnote{We use a different naming convention for our Ce-vacancy complexes because the Ce atoms do not remain at substitutional sites after relaxation, as implied by the naming convention of Ref.~\cite{Tan2020Ce}: CeV, CeV$2$ and CeV$3$.}

We performed DFT calculations in \texttt{CASTEP} \cite{Clark2005} on the input structures, which were randomly distorted from ideal symmetry to ensure energetically most favorable final structures after optimization, using on-the-fly-generated (OTFG) ultra-soft pseudopotentials \cite{vanderbilt1990}. An energy cutoff of $540$ eV and a $3 \times 3 \times 3$ reciprocal-space sampling k-mesh were employed. To account for localisation of the $4f$ orbitals on Ce, we used PBE$+U$ \cite{Perdew1996,Dudarev1998} with the default value of $U=6.0$ eV for Ce in CASTEP. The same range of initial charge states as in Ref.~\cite{Tan2020Ce} were used: $q=-1,0,+1,+2,+3,+4$. In addition, to be able to account for all possible defect spin states, different initial spin states for each charge state were considered: $0,1,2,...,6,7$. The structures were relaxed to tolerance $0.01$ eV/\AA. \footnote{The \texttt{CASTEP} utility \texttt{c2x} was employed to convert output binary files to Gaussian \texttt{cube} format read by \texttt{SPINNEY}.}

\begin{figure}[htb]
\centering
\includegraphics[width=0.3095\textwidth]{DIA215-CeV-Rand_initial}
\includegraphics[width=0.319\textwidth]{DIA214-CeV2-Rand_initial}
\includegraphics[width=0.3205\textwidth]{DIA213-CeV3-Rand_initial}
\caption{\label{fig1} The initial (\emph{left}) Ce$_{\rm V2}$, (\emph{middle}) Ce$_{\rm V3}$ and (\emph{right}) Ce$_{\rm V4}$ complexes used in this work. Faded yellow stars indicate the nearest-neighbour atoms of Ce that have been removed by analogy with Ref.~\cite{Tan2020Ce}.}
\end{figure}

% \begin{figure}[htb]
% \centering
% \includegraphics[width=0.3575\textwidth]{DIA215-CeV-Rand_opt_atoms}
% \includegraphics[width=0.3035\textwidth]{DIA214-CeV2-Rand_opt_atoms}
% \includegraphics[width=0.2985\textwidth]{DIA213-CeV3-Rand_opt_atoms}
% \caption{\label{fig2} (The final (\emph{left}) Ce$_{\rm V2}$, (\emph{middle}) Ce$_{\rm V3}$ and (\emph{right}) Ce$_{\rm V4}$ vacancy complexes. Only the dopant and nearest-neighbour carbon atoms within the second nearest-neighbour shelll of pristine diamond are shown. The bond distance and most angles are also shown.}
% \end{figure}

\begin{figure}[htb]
\centering
\includegraphics[width=0.31\textwidth]{DIA215-CeV-Rand_opt}
\includegraphics[width=0.32\textwidth]{DIA214-CeV2-Rand_opt}
\includegraphics[width=0.33\textwidth]{DIA213-CeV3-Rand_opt}
\caption{\label{fig2} The optimized (\emph{left}) Ce$_{\rm V2}$, (\emph{middle}) Ce$_{\rm V3}$ and (\emph{right}) Ce$_{\rm V4}$ complexes. The dopant and nearest-neighbour ($nn$) carbon atoms within the second $nn$ shell of pristine diamond are shown. The bond distances are also shown. In Ce$_{\rm V2}$, all six nearest neighbours are $\approx2.2\,\text{\AA}$ from the Ce atom. In Ce$_{\rm V3}$, two nearest neighbours are $\approx2.2\,\text{\AA}$, and the other four $\approx2.3\,\text{\AA}$, from the Ce atom. In Ce$_{\rm V4}$, one nearest neighbour is $\approx2.2\,\text{\AA}$, and the other three $\approx2.3\,\text{\AA}$, from the Ce atom. All other carbon atoms are further than the second $nn$ shell of pristine diamond ($2.53\,\text{\AA}$) from the Ce atom.}
\end{figure}

% Try to avoid unnecessary white space as much as possible.  When making a new paragraph within the same section or subsection, you should add a blank line like this to space new paragraphs apart, as seen here.\\ %skip a line when changing paragraph, the double \\ and a few line changes in the code like this.

% {\bf Small Bug:} if your affiliations are squashed at the top below the authors' names on the first page when compiling (shown below), try to use a \textbackslash newpage to force a page change earlier/later to avoid the problem.  

% \includegraphics[width=0.999\textwidth]{squashed}

% put a \newpage on page 1 before/after the automatic page change to the second page if your affiliations are squashed.  If you don't understand what this is all about, and your document looks fine, just ignore this statement.
%\newpage

\section{Results}
The relaxed vacancy complexes are shown in Fig.~\ref{fig2}. The dopant and nearest neighbour carbon atoms within the second-nearest neighbour shell of pristine diamond  ($2.53\,\text{\AA}$) are highlighted. The extent of coordination and the increasing disparity between neighbouring bond lengths indicate the growing disorder about the relaxed point defects.

It should not come as a surprise that the FNV corrections become increasingly unreliable as the disorder about the dopant increases \cite{Kumagai2014}. This can be seen by comparing the obtained overall FNV and KO corrections (in eV) in tables~\ref{table1}--\ref{table3}. We also calculated the enclosed charge within each supercell based on the fitting parameters $x$, $\beta$ and $\gamma$ obtained from Eq.~(\ref{eq:FNVrho}) during the FNV correction calculation. Then, to enable a comparison with the conventional potential alignment term in the second-to-last column, the FNV and KO potential alignment terms are also listed. To obtain the conventional values, we compared the electrostatic potential of the defect supercells at an interstitial centroid far away from the defect with that of the relaxed pristine supercell at the same location. 

\begin{table}[htb!]
\centering % centering table
\setlength{\tabcolsep}{0.5em} % for the horizontal padding
%{\renewcommand{\arraystretch}{1.2}% for the vertical padding
\begin{tabular}{r c r r r r r r}
\hline\hline
$q$ & $S_{\min}$ &
\multicolumn{1}{c}{FNV total} &
\multicolumn{1}{c}{FNV $\Delta V_{\text{far}}$} &
\multicolumn{1}{c}{KO total} &
\multicolumn{1}{c}{KO $\Delta V_{\text{far}}$} &
\multicolumn{1}{c}{Conventional $\Delta V$} &
 \multicolumn{1}{c}{ $q_{\mathrm{enc}}$} \\ & & \multicolumn{5}{c}{(eV)}\\
\hline
$-1$ & $0.98754$ & $0.226$ & $0.118$ & $0.217$ & $0.126$ & $0.163$ & $-0.99999996$\\
$0$  & $0.0000022$ & $0.0$ & $-0.003$ & $0.0$ & $0.001$ & $0.0$ & $0.0$\\
$+1$ & $1.0007$ & $0.205$ & $-0.119$ & $0.210$ & $-0.114$ & $-0.157$ & $0.99999994$\\
$+2$ & $2.00004$ & $0.774$ & $-0.222$ & $0.771$ & $-0.222$ & $-0.303$ & $1.9999998$\\
$+3$ & $2.08479$ & $1.893$ & $-0.346$ & $1.796$ & $-0.345$ & $-0.3134$ & $2.9999994$\\
$+4$ & $2.17017$ & $27.245$ & $-0.098$ & $0.865$ & $-0.087$ & $-0.3365$ & $3.9998$ \\
\hline
\end{tabular}%}
\caption{\label{table1}FNV and KO total correction and potential alignment values (eV) for vacancy complex Ce$_{\rm V2}$. $S_{\min}$ corresponds to the spin state of lowest total energy. The conventional potential alignment is included for comparison. The enclosed charge within each supercell is calculated as $q_{\mathrm{enc}}(r)
  = 4\pi \int_{0}^{r} \rho_{\rm model,d}(r')\, r'^{2}\,dr'$ with $\rho_{\rm model,d}$ as defined in Eq.~(\ref{eq:FNVrho}) and the integration performed on interval 
$0 < r \le \tfrac12 L_{\min}$ with $L_{\min} = \min_{i}\,\lVert \mathbf a_{i}\rVert$ and $\mathbf a_{i}$ the supercell side lengths.} 
\end{table}

\begin{table}[htb!]
\centering % centering table
\setlength{\tabcolsep}{0.5em} % for the horizontal padding
%{\renewcommand{\arraystretch}{1.2}% for the vertical padding
\begin{tabular}{r c r r r r r r}
\hline\hline
$q$ & $S_{\min}$ &
\multicolumn{1}{c}{FNV total} &
\multicolumn{1}{c}{FNV $\Delta V_{\text{far}}$} &
\multicolumn{1}{c}{KO total} &
\multicolumn{1}{c}{KO $\Delta V_{\text{far}}$} &
\multicolumn{1}{c}{Conventional $\Delta V$} &
 \multicolumn{1}{c}{ $q_{\rm enc}$} \\ & & \multicolumn{5}{c}{(eV)}\\
\hline
$-1$ & $1.0$ & $0.317$ & $0.113$ & $0.202$ & $0.116$ & $0.121$ & $-0.9999995$\\
0 & $0.0019$ & $0.0$ & $-0.003$ & $0.0$ & $0.003$ & $0.0$ & $0.0$\\[2pt]
$+1$ & $0.99811$ & $297.953$ & $297.953$ & $0.201$ & $-0.102$ & $-0.134$ & $7.3 \times 10^{-6}$ \\
$+2$ & $0.00062$ & $1498.712$ & $749.356$ & $0.725$ & $-0.197$ & $-0.227$ & $9.2\times 10^{-6}$ \\
$+3$ & $0.787866$ & $2967.798$ & $989.266$ & $1.415$ & $-0.268$ & $-0.288$ & $1.4\times 10^{-5}$ \\
$+4$ & $0.188895$ & $3526.123$ & $881.531$ & $1.712$ & $-0.234$ & $-0.345$ & $1.8\times 10^{-5}$ \\
\hline
\end{tabular}%}
\caption{\label{table2} FNV and KO total correction and potential alignment values (eV) for vacancy complex Ce$_{\rm V3}$. The meanings of the remaining columns carry over from table~\ref{table1}.} 
\end{table}

\begin{table}[htb!]
\centering % centering table
\setlength{\tabcolsep}{0.5em} % for the horizontal padding
% {\renewcommand{\arraystretch}{1.2}% for the vertical padding
\begin{tabular}{r c r r r r r r}
\hline\hline
$q$ & $S_{\min}$ &
\multicolumn{1}{c}{FNV total} &
\multicolumn{1}{c}{FNV $\Delta V_{\text{far}}$} &
\multicolumn{1}{c}{KO total} &
\multicolumn{1}{c}{KO $\Delta V_{\text{far}}$} &
\multicolumn{1}{c}{Conventional $\Delta V$} &
 \multicolumn{1}{c}{ $q_{\rm enc}$} \\ & & \multicolumn{5}{c}{(eV)}\\
\hline
$-1$ & $1.0$ & $292.236$ & $-292.236$ & $0.147$ & $0.079$ & $0.119$ & $-4.6\times 10^{-6}$\\
$ 0$ & $0.0$ &   $0.0$ & $-0.001$   & $0.0$ & $0.001$ & $0.0$ & $0.0$   \\
$+1$ & $0.990084$ & $285.872$ & $285.872$ & $0.186$ & $-0.104$ & $-0.122$ & $7.3\times 10^{-6}$\\
$+2$ & $0.00087$ & $1113.150$ & $556.575$ & $0.727$ & $-0.202$ & $-0.237$ & $1.4\times 10^{-5}$ \\
$+3$ & $0.961946$ & $2218.033$ & $739.344$ & $1.450$ & $-0.268$ & $-0.322$ & $2.2\times 10^{-5}$\\
$+4$ & $0.337866$ & $2843.334$ & $710.833$ & $1.886$ & $-0.252$ & $-0.392$ & $2.9\times 10^{-5}$\\
\hline
\end{tabular}%}
\caption{\label{table3}FNV and KO total correction and potential alignment values (eV) for vacancy complex Ce$_{\rm V4}$. The meanings of the remaining columns carry over from table~\ref{table1}.} 
\end{table}

With the exception of charge state $q=+4$, Ce$_{\rm V2}$ in table~\ref{table1} exhibits consistent values of the total FNV and KO corrections. The potential alignments of the two schemes also agree and both agree reasonably well with the conventional scheme. For the $q=+3$ charge state, possible metastable spin states have been found; only the most stable spin states are reported: $S_{\min}$ corresponds to the spin state of lowest total energy. The other spin-state total energies of $q=+3$ differ from the minimum by at most $6.0\times 10^{-5}$ eV. Hence, given this small value, we conclude that there are no metastable spin states for this vacancy complex. Finally, the FNV model charge calculated from Eq.~(\ref{eq:FNVrho}) for each defect supercell agrees with the charge state specified in the CASTEP calculation. The discrepancy between the FNV and KO corrections for charge state $q=+4$ suggests that it is highly anisotropic. This is an example of the inherit deficiency in the FNV method that the KO correction was meant to overcome\cite{Kumagai2014}. 

We immediately observe how the situation changes when we consider the next vacancy complex with increased disorder: Ce$_{\rm V3}$. FNV only yields reasonable total correction and potential alignment values for the $q=-1$ charge state. The FNV model charge enclosed in the supercell also tallies with this charge. The FNV correction and potential alignment values blow up for all positive charge states, and the fact that the essentially zero fitted charge densities do not agree with them suggests that the charge is delocalised, and that these defects are shallow, unlike the localised deep states in Ce$_{\rm V2}$. For $q=0$, all non-zero initial spin states turn out to be metastable, with total energies that differ by up to $\approx 0.01$ eV from $S_{\min}$. These energetically less favourable spin states have $S\approx2.0$ compared to $S_{\min}\approx0.002$. 

For the largest and most disordered vacancy cluster, Ce$_{\rm V4}$, the FNV correction fails across the board, and all non-zero charge states yield shallow defects with the FNV model charge delocalised beyond the confines of a single defect supercell. The total energy difference between metastable spin states (for a given charge state) is up to $\approx 0.05$ eV for $q=+2$. The less energetically favourable spin states all end up as $S\approx2.0$ whilst $S_{\min}\approx1\times 10^{-5}$. In constrast to $q=+4$ of Ce$_{\rm V2}$, the FNV correction ultimately fails here not because the defect is anisotropic but rather because it is shallow, as is evident from the delocalised model charge in the last columns of tables~\ref{table2} and \ref{table3}.

The KO correction, on the other hand, delivers consistent values across all charge states. This is to be expected given that it employs a point-charge electrostatic and atomic-site potential alignment model which by design cannot discern the spatial extent of the charge distribution. Where KO agrees with the FNV correction, we know that we are dealing with deep defect states

Finally, it is worth mentioning that the stability of the neutral vacancy complexes in our work follows the same trend as the relaxed defects in Refs.~\cite{Tan2020Ce,Morris2024}: $E_{f}\!\bigl[D^{0}\bigr]_{\rm Ce_{\rm V2}}=15.5259$ eV, $E_{f}\!\bigl[D^{0}\bigr]_{\rm Ce_{\rm V3}}=14.3588$ eV and $E_{f}\!\bigl[D^{0}\bigr]_{\rm Ce_{\rm V4}}=14.5393$ eV. In other words, the neutral $\rm Ce_{\rm V3}$ complex is the most stable. 

In conclusion, the FNV and KO correction schemes have been evaluated for a set of three Ce-vacancy complexes across various charge states, and spin states for each charge. The FNV correction fails gracefully in the presence of a shallow defect whilst the KO correction continues to provide reasonable values. FNV results in unrealistically large correction values for shallow levels, allowing defects with shallow levels to be identified.  Ultimately, the cure for a failing FNV correction is to employ hybrid exchange-correlation functionals as is our intention in future work. 

%\clearpage

%Bibliography
%make sure you fill the .bib properly, e.g. edit "biblio.bib" and add all your references in there.
\bibliographystyle{IEEEtran}
\bibliography{biblio}

\end{document}