\section{Introduction}
Two-dimensional (2D) materials provide the unique opportunity to scale future electronics smaller than ever believed physically possible, implying engineering 2D materials is a promising strategy that can meet the demands of future nanotechnologies \cite{Butler2013}.
As defects play a crucial role in the optical and electronic properties of these systems, engineering of defects in 2D materials has sparked continuous interest \cite{wang,hong2017atomic,Lin2016,Dreyer2018,weston2018native,Tawfik2017}. For example, defects in \textit{h}-BN have been found to be the source of stable polarized and ultra-bright single-photon emissions at room temperature\cite{bourrellier2016bright,Abdi2018,Aharonovich170,Tran2015}. Hence, the development of our understanding of defects in 2D materials will open up further possibilities for emerging applications in quantum information and nanotechnology with much better scalability than traditional defects in 3D materials.

Unlike in their 3D counterparts \cite{Freysoldt2009,Freysoldt2014,Vinichenko2017,Kumagai2014,Komsa2013}, first-principles techniques for calculating defect properties in 2D materials still face significant challenges.  Specifically, eliminating the periodic charge interactions for charged defects in 2D materials requires a charge correction scheme that accounts for the weak and anisotropic dielectric screening of 2D systems \cite{Komsa2014,Wang2015}. Furthermore, several exchange-correlation functionals that provide accurate electronic structures for 3D bulk systems are no longer applicable to ultrathin 2D systems.  
For example, the fraction of Fock exchange ($\alpha$) in hybrid functionals can be approximated as the inverse of dielectric constant ($\varepsilon_{\infty}$) of the material \cite{alkauskas2011defect,Jonathan2014}; however, $\varepsilon_{\infty}$ cannot be uniquely defined for low dimensional systems as discussed in previous studies.\cite{Delerue2003,Hamel2008,brawand2016generalization}
Therefore, the determination of $\alpha$ in hybrid functionals for 2D materials remains an open question. On the other hand, many body perturbation theory techniques ($\textit{e.g.}$ GW approximation) give accurate quasiparticle energies such as band gaps and band positions; however, high computational cost and slow convergence with respect to empty states make the screening of many defects in 2D materials impractical with conventional implementations \cite{Qiu2016,Rasmussen2016,Thygesen2017,Huser2013,Attaccalite2011,Jornada2017}.

In our previous work \cite{feng,Sundararaman2017}, we developed an efficient and accurate method that can give reliable charge corrections for total energies and electronic states of charged defects in 2D materials \textit{without any supercell extrapolations}, and then provided accurate defect CTLs with the DFT+GW scheme \cite{Malashevich2014,Chen2013,Chen2015,Chen2017}.
Such  implementation is built on top of the WEST-code \cite{west}, Quantum-Espresso \cite{Giannozzi2017} and JDFTx \cite{Sundararaman2017b} packages. In our GW calculations, we avoided explicit inclusion of empty states and inversion of dielectric matrices \cite{west, ping2013,anh2013}, while also speeding up vacuum size convergence with a 2D Coulomb truncation \cite{Ismail-Beigi2006}. 

In this paper we will further investigate the possibilities of even less expensive but similarly reliable computational methods for both defect charge transition levels and ionization energies in 2D materials. The rest of the paper is organized as follows. In section II, we discuss the methodology, including details of computational methods in II.A. and how to compute thermodynamic charge transition levels in general in II.B. Then in section III we discuss our results where we have addressed two important issues for 2D materials: in section III.A., we determine which level of theory and which electron chemical potential reference should be used to calculate a CTL in 2D systems; in section III.B., we show how to define the fraction of Fock exchange in hybrid functionals for accurate band edges and band gaps; in the end, section III.C., we combine these two findings to obtain accurate defect ionization energies for 2D materials.


\section{Methodology}

\subsection{Computational Methods}
In this work, all structural relaxations and total energy calculations were performed using open source plane wave code Quantum-ESPRESSO \cite{Giannozzi2009} with ONCV norm-conserving pseudopotentials \cite{Hamann2013,Schlipf2015}, a wavefunction cutoff of 70 Ry and a $k$-point mesh equivalent to $12\times12\times1$ or higher in the primitive cell. 
 The GW calculations were performed using the WEST code \cite{Govoni2015}, which avoids explicit empty states and inversion of dielectric matrices. We employed the contour deformation technique for frequency integration of the self energy. The final values of GW corrections were extrapolated between $9\times 9$ and $12\times 12$ $k$-point meshes to infinite $k$-points similar to Ref.~\citenum{feng}. A two-dimensional Coulomb truncation~\cite{Ismail-Beigi2006,feng} has been applied to speed up the vacuum size convergences. The charge corrections for the total energies and eigenvalues of charged defects employed the techniques developed in Ref. \citenum{Sundararaman2017,feng}, which were implemented in the JDFTx code \cite{Sundararaman2017b,Ismail-Beigi2000,Arias1992}. More computational details can be found in supplementary materials.

\subsection{Thermodynamic Charge Transition Levels}

\begin{figure}[t]
\includegraphics[keepaspectratio=true,width=\linewidth]{path.png}
\caption{Schematic plot of the two paths (distinguished with blue/red color) that transition from charge state $q$ to $q+1$. For each path, there is a corresponding vertical excitation, which can be computed either with EA$_{q+1}$ or IP$_q$ (noted with up/down arrowheads), as discussed in the main text.}  \label{fig:path}
\end{figure}


A thermodynamic CTL is the value of electron chemical potential $\varepsilon_F$ at which the stable charge state of the system changes, \textit{e.g.} from $q$ to $q+1$. Therefore, CTLs are calculated through the equivalency of the formation energies $q$ and $q+1$, given by Eq. (\ref{eq:ctl0}) \cite{Freysoldt2014}.
\begin{align}
\varepsilon_{q+1|q}&=E_q^f(\mathbf{R}_q)-E_{q+1}^f(\mathbf{R}_{q+1}) \nonumber\\
&=E_q(\mathbf{R}_q)-E_{q+1}(\mathbf{R}_{q+1})-\varepsilon_F
\label{eq:ctl0}
\end{align}
Here $E_q^f(\mathbf{R})$ is the defect formation energy with charge $q$ and geometry $\mathbf{R}$, and $\mathbf{R}_q$ is the relaxed geometry of the system with charge $q$. $E_q(\mathbf{R})$ is the total energy that relates to $E_q^f(\mathbf{R})$ and $\varepsilon_F$ following the definition of Eq. (1) in Ref. \citenum{feng}.  Diagrammatically, Eq. (\ref{eq:ctl0}) is the energy difference between two potential  surface minima in position space $\mathbf{R}$, as shown in Fig. \ref{fig:path}.  



\section{Results and Discussions}
\subsection{Implementing Quasiparticle Corrections in Defect Charge Transition Levels}

In extended systems, local and semi-local functionals fail to yield accurate total energy differences between two charge states, i.e. where an electron removing (IP)/adding process (EA) is involved. An alternative approach \cite{feng} is to separate Eq.~(\ref{eq:ctl0}) into two parts: the vertical excitation energy between two charge states ($q$ and $q+1$) at the same geometry ($\mathbf{R}$) (denoted as quasiparticle energies $\varepsilon^{QP}$) and the geometry relaxation energy at a fixed charge state (denoted with $E^{rlx}$). 
In general, as ground states theory, DFT provides reliable geometry relaxation energies at a fixed charge state in many systems. For the 2D systems we tested in this work, we found the total energy difference of geometries optimized at semi-local and hybrid functionals is less than 10 meV, which leads to a negligible difference in charge transition levels (see SI Table IV). Hence, this separation allows us to accurately calculate the vertical excitation energies with a higher level of theory appropriate for non-neutral excitations, such as the GW approximation.

One can separate Eq.~(\ref{eq:ctl0}) by two possible physical pathways from $E_q^f(\mathbf{R}_q)$ to $E_{q+1}^f(\mathbf{R}_{q+1})$ as shown in Fig. \ref{fig:path}. One pathway (red path) occurs with a vertical excitation at $\mathbf{R}_{q}$ ($E_{q+1}^f(\mathbf{R}_q)-E_q^f(\mathbf{R}_q)$) followed by a geometry relaxation at the charge state $q+1$ ($E_{q+1}^f(\mathbf{R}_{q+1})-E_{q+1}^f(\mathbf{R}_q)$), shown in Eq.~(\ref{eq:ctl1}). The other pathway (blue path) occurs through the geometry relaxation at the charge state $q$ plus a vertical excitation at $\mathbf{R}_{q+1}$, corresponding to Eq.~(\ref{eq:ctl2}).
\begin{align}
\varepsilon_{q+1|q} &= \underbrace{E_{q}^f(\mathbf{R}_q)-E_{q+1}^f(\mathbf{R}_q)}_{\varepsilon^{QP}}+\underbrace{E_{q+1}^f(\mathbf{R}_q)-E_{q+1}^f(\mathbf{R}_{q+1})}_{E^{rlx}}  \nonumber\\
&= \varepsilon^{QP}_{q+1|q}(\mathbf{R}_q) + E^{rlx}_{q+1} 
\label{eq:ctl1}
\end{align}
\begin{align}
\varepsilon_{q+1|q}&=\underbrace{E_q^f(\mathbf{R}_q)-E_q^f(\mathbf{R}_{q+1})}_{E^{rlx}}+\underbrace{E_q^f(\mathbf{R}_{q+1})-E_{q+1}^f(\mathbf{R}_{q+1})}_{\varepsilon^{QP}} \nonumber\\
&= E^{rlx}_q  + \varepsilon^{QP}_{q+1|q}(\mathbf{R}_{q+1})
\label{eq:ctl2}
\end{align}
Note that all three equations (Eq. (\ref{eq:ctl0}), (\ref{eq:ctl1}), (\ref{eq:ctl2})) are theoretically equivalent. Yet, in practice they may yield sizable differences, when $\varepsilon^{QP}$ is computed through eigenvalues at a specific level of theory instead of total energy differences.

Furthermore, the vertical excitation energies $\varepsilon^{QP}_{q+1|q}$ in Eq.~(\ref{eq:ctl1}) and Eq.~(\ref{eq:ctl2}) can be determined from either the ionization potential of the charge state $q$ ($\text{IP}_q$) or
the electron affinity of the charge state $q+1$ ($\text{EA}_{q+1}$)
, as noted in Fig.~\ref{fig:path} with up/down arrowheads. We obtained IP and EA through eigenvalues at different levels of theory based on the Janak's theorem \cite{Janak}.
Within the framework of DFT with local and semi-local functionals, energy is a continuous and smooth functional of the number of electrons in the system $E[n]$. This results in non-linear behavior with respect to electron number and gives inconsistent eigenvalues ($\partial E / \partial n$) of the $q$ and $q+1$ systems (resulting in $\text{IP}\neq\text{EA}$) \cite{perdew1982density,cohen2008insights}. Following the discussions in Ref. \citenum{Fabien2009}, an eigenvalue between the $q$ and the $q+1$ systems can be approximated to the second order as:
\begin{equation}
  \varepsilon^{QP}_{q+1|q}(\mathbf{R}) \approx - \frac{1}{2}\left[ \frac{\partial E}{\partial n} \bigg|_{n=(q+1)^-} + \frac{\partial E}{\partial n} \bigg|_{n=q^+} \right]
  \label{eq:avg1}
\end{equation}
Therefore, to the second order in electron number, $\varepsilon^{QP}$ is obtained by taking the average of $\text{IP}_q$ and $\text{EA}_{q+1}$:
\begin{equation}
  \varepsilon^{QP}_{q+1|q}(\mathbf{R}) = \frac{1}{2}\bigg[\text{EA}_{q+1}(\mathbf{R}) + \text{IP}_q(\mathbf{R})\bigg]
  \label{eq:avg2}
\end{equation}


\begin{table}[t]
\centering
\comment{\begin{tabular}{p{1.5cm}p{1.5cm}p{1.2cm}p{1.2cm}p{1.2cm}p{1.2cm}}
\hline\hline
\multicolumn{2}{c}{Method} & \multicolumn{4}{c}{Defect} \\
 & & C$_\text{B}$  & V$_\text{N}$C$_\text{B}$  & C$_\text{N}$ & V$_\text{N}$C$_\text{B}$  \\
\hline
&&&&&\\
     &    CTL             & (0/+1) & (0/+1) & (-1/0) & (-1/0) \\
   &&&&&\\   
     & Eq\ref{eq:ctl0} & -3.63  & -4.22  & -3.54  & -1.57  \\
PBE  & Eq\ref{eq:ctl1} & -3.61  & -4.29  & -3.51  & -1.66  \\
     & Eq\ref{eq:ctl2} & -3.64  & -4.33  & -3.49  & -1.67  \\
   
 &&&&&\\
     & Eq\ref{eq:ctl0} & -3.65  & -4.19  & -3.50  & -1.87  \\
PBE0 & Eq\ref{eq:ctl1} & -3.60  & -4.17  & -3.50  & -1.87  \\
     & Eq\ref{eq:ctl2} & -3.62  & -4.21  & -3.50  & -1.21*  \\
   
&&&&&\\
     & Eq\ref{eq:ctl1} & -3.40  & -4.29  & -3.74  & -1.74 \\
$G_0W_0$   & Eq\ref{eq:ctl2} & -3.28  & -4.22  & -3.73  & -1.70  \\
   
    &&&&&\\
    \multicolumn{6}{c}{IP$_{q}(\mathbf{R}_{q})$-EA$_{q+1}(\mathbf{R}_{q})$} \\
     &&&&&\\
PBE && 2.68  & 2.60  & 2.75  & 2.50  \\
PBE0 && 1.15  & 1.09  & 1.13  & 1.42  \\
$G_0W_0$  && 0.04  & 0.20  & 0.03  & 0.19  \\
 \hline\hline
\end{tabular}}
\begin{tabular}{p{1cm}p{1.1cm}p{1.1cm}p{1.1cm}p{1.1cm}p{1.1cm}| @{\hspace{2mm}} c}
\hline\hline
\multicolumn{2}{c}{ } & \multicolumn{5}{c}{Defect} \\
\multicolumn{2}{c}{Method} & \hspace{2mm}C$_\text{B}$  & V$_\text{N}$C$_\text{B}$  & \hspace{2mm}C$_\text{N}$ & V$_\text{N}$C$_\text{B}$  & V$_\text{S}$ \\
\multicolumn{2}{c}{  } &&\multicolumn{2}{c}{(\textit{h}-BN)} && (MoS$_2)$\\
\hline
&&&&&\\
     &    CTL             & (0/+1) & (0/+1) & (-1/0) & (-1/0) & (-1/0) \\
   &&&&&\\  
PBE  & Eq\ref{eq:ctl0} & -3.63  & -4.22  & -3.54  & -1.57 & -4.29  \\
     & Eq\ref{eq:ctl1} & -3.61  & -4.29  & -3.51  & -1.66 & -4.29  \\
     & Eq\ref{eq:ctl2} & -3.64  & -4.33  & -3.49  & -1.67 & -4.29  \\
   
 &&&&&\\
PBE0 & Eq\ref{eq:ctl0} & -3.65  & -4.19  & -3.50  & -1.87 & -4.33  \\
     & Eq\ref{eq:ctl1} & -3.60  & -4.17  & -3.50  & -1.87 & -4.32  \\
     & Eq\ref{eq:ctl2} & -3.62  & -4.21  & -3.50  & -1.21* & -4.31  \\
   
&&&&&\\
$G_0W_0$   & Eq\ref{eq:ctl1} & -3.40  & -4.29  & -3.74  & -1.74 & -4.34 \\
     & Eq\ref{eq:ctl2} & -3.28  & -4.22  & -3.73  & -1.70 & -4.38  \\
   
    \multicolumn{7}{c}{ }\\
    \multicolumn{7}{c}{IP$_{q}(\mathbf{R}_{q})$-EA$_{q+1}(\mathbf{R}_{q})$} \\
    \multicolumn{7}{c}{ }\\
PBE && 2.68  & 2.60  & 2.75  & 2.50 & 0.947  \\
PBE0 && 1.15  & 1.09  & 1.13  & 1.42 & 0.023  \\
$G_0W_0$  && 0.04  & 0.20  & 0.03  & 0.19 & -0.202  \\
 \hline\hline
\end{tabular}
\caption{\label{table:path} Charge transition levels (CTLs) relative to vacuum (in eV) of multiple defects in monolayer \textit{h}-BN and V$_\text{S}$ in monolayer MoS$_2$. These values are collected via three methods (Eq. (\ref{eq:ctl0}-\ref{eq:ctl2})) at various levels of theory (PBE, PBE0, $G_0W_0$$@$PBE ). The CTLs relative to vacuum are remarkably similar. The one exception, V$_\text{N}$C$_\text{B}$  (-1/0) at PBE0 (marked with *) incidentally has a band inversion resulting in a CTL within the valence band, breaking the reliability of Eq. (\ref{eq:ctl2}).
We also show IP$_{q}(\mathbf{R}_{q})-$EA$_{q+1}(\mathbf{R}_{q})$ at different levels of theory. Note that at the $G_0W_0$  level, this difference is $<0.2$ eV.}
\end{table} 


By employing this principle in Eq.~\ref{eq:avg2}, we compared the CTL obtained with PBE, PBE0 and $G_0W_0@$PBE for three different defects in monolayer BN
and S vacancy in MoS$_2$ as shown in Table \ref{table:path}.
Here we propose to set $\varepsilon_F$ equal to the vacuum level (determined by the electrostatic potential in the vacuum region of supercells) and use it as a reference for Eq.(\ref{eq:ctl0}).
We found this choice (opposed to the commonly used band edges) is particularly advantageous for obtaining consistent CTLs among different methods as shown in Table~\ref{table:path}.
(More computational details for $G_{0}W_{0}$ can be found in SI, with numerical techniques as in Ref. \citenum{feng}).
We note that our choice of using the vacuum level as the reference for 2D materials has similarity with the idea of using the averaged electrostatic potential as the reference level for 3D materials, as discussed in several previous studies.~\cite{alkauskas2008defect,lyons2009role,komsa2011assessing,peng2013convergence}
There are several interesting observations from Table \ref{table:path}, as follows. First, we found excellent agreement (within 0.1 eV) among Eq.~(\ref{eq:ctl0}), (\ref{eq:ctl1}) and (\ref{eq:ctl2}) for each defect at a fixed level of theory if $\varepsilon^{QP}_{q+1|q}$ in Eq.~(\ref{eq:ctl1}) and (\ref{eq:ctl2}) is obtained through Eq.~(\ref{eq:avg2}). Second, we found the results obtained among PBE, PBE0 and $G_0W_0@$PBE are also strikingly similar (within 0.2 eV) for several defects in \textit{h}-BN as well as 
V$_\text{S}$ in MoS$_2$ which has a very different chemical bonding character from \textit{h}-BN.
This suggests that CTLs of 2D materials relative to vacuum are \textit{robust to the level of theory one chooses}. Note that the difference between $\text{IP}_q$ and $\text{EA}_{q+1}$ is more than 2 eV for PBE, reduced to $~$1 eV at PBE0 level ($\alpha=0.25$), but less than 0.2 eV at $G_0W_0$, which indicates the delocalization error present in semi-local DFT has been mostly corrected at $G_0W_0$@PBE \cite{Fabien2009}. To the best of our knowledge, there are currently no available experimental values of defect charge transition levels in ultrathin 2D materials to compare with our calculations, partially due to the difficulty of controlling and identifying the chemical composition of defects in 2D materials. However, this emphasizes the need to reliably predict defect charge transition levels in 2D materials 
from first principles, and then use them to identify the chemical composition of defects by comparing with experimental measurements, such as zero phonon lines (ZPL). We note that we compared our results with the previous theoretical studies at the corresponding level of theory referenced to the valence band maximum (VBM), and obtained overall good agreement for monolayer and bulk \textit{h}-BN as well as MoS$_2$, as shown in SI Table V. 

\begin{table*}[t]
\begin{center}
\begin{tabular}{p{1.9cm}p{2.1cm}p{2.1cm}p{2.1cm}p{2.1cm}p{2.1cm}p{2.1cm}p{2cm}}
	 \hline\hline
	 \hspace{2mm} System & \hspace{4mm} PBE & \hspace{4mm} HSE & \hspace{3mm} PBE0 & \hspace{2mm} B3PW & \hspace{1mm} PBE0($\alpha$)\text{} & \hspace{3mm} $G_0W_0$  &
	 \hspace{2mm} EXP.\\
	 \hline
	\hspace{0.1mm} ML MoS$_2$
	& 1.74 $|$ K$\,\rightarrow\,$K & 2.17 $|$ K$\,\rightarrow\,$K
	& 2.85 $|$ K$\,\rightarrow\,$K & 2.58 $|$ K$\,\rightarrow\,$K
	& 2.85 $|$ K$\,\rightarrow\,$K & 2.82 $|$ K$\,\rightarrow\,$K & 2.7\cite{Krane2016MoS2,Yao2017MoS2} \\
	\hspace{0.1mm} Graphane
	& 3.57 $|$ $\Gamma$$\,\rightarrow\,$$\Gamma$  & 4.41 $|$ $\Gamma$$\,\rightarrow\,$$\Gamma$ 
	& 5.06 $|$ $\Gamma$$\,\rightarrow\,$$\Gamma$  & 5.04 $|$ $\Gamma$$\,\rightarrow\,$$\Gamma$ 
	& 6.54 $|$ $\Gamma$$\,\rightarrow\,$$\Gamma$  & 6.41 $|$ $\Gamma$$\,\rightarrow\,$$\Gamma$ 
	& -- \\
	\hspace{0.1mm} ML BN 
	& 4.71 $|$ K$\,\rightarrow\,$K & 5.70 $|$ K$\,\rightarrow\,$$\Gamma$ 
	& 6.33 $|$ K$\,\rightarrow\,$$\Gamma$  & 6.33 $|$ K$\,\rightarrow\,$$\Gamma$ 
	& 7.34 $|$ K$\,\rightarrow\,$$\Gamma$  & 7.01 $|$ K$\,\rightarrow\,$$\Gamma$   & -- \\ 
	\hspace{0.1mm} BL BN 
	& 4.49 $|$ T$\,\rightarrow\,$M & 5.81 $|$ T$\,\rightarrow\,$M
	& 6.46 $|$ T$\,\rightarrow\,$$\Gamma$  & 6.17 $|$ T$\,\rightarrow\,$M
	& 7.08 $|$ T$\,\rightarrow\,$$\Gamma$  & 7.00 $|$ T$\,\rightarrow\,$$\Gamma$   & -- \\ 
	\hspace{0.1mm} TL BN 
	& 4.36 $|$ T$\,\rightarrow\,$M & 5.68 $|$ T$\,\rightarrow\,$M
	& 6.40 $|$ T$\,\rightarrow\,$M & 6.03 $|$ T$\,\rightarrow\,$M
	& 7.01 $|$ T$\,\rightarrow\,$$\Gamma$  & 6.92 $|$ T$\,\rightarrow\,$M & -- \\ 
	\hspace{0.1mm} Bulk BN
	& 4.22 $|$ T$\,\rightarrow\,$M & 5.60 $|$ T$\,\rightarrow\,$M
	& 6.28 $|$ T$\,\rightarrow\,$M & 5.91 $|$ T$\,\rightarrow\,$M
	& 6.07 $|$ T$\,\rightarrow\,$M & 6.01 $|$ T$\,\rightarrow\,$M  &6.08 $\pm$ 0.015 \\	
	 \hline
	 \hline
	\end{tabular}
      \captionof{table}{Electronic band gaps (eV) for various pristine 2D materials. In general, PBE severely underestimates the gap. Hybrid functionals HSE, B3PW, and PBE0 ($\alpha=0.25$) generally enlarge the bulk band gap, but still underestimate the gaps of ultrathin BN and graphane compared with experiments and GW approximation. Only PBE0($\alpha$) with $\alpha$ satisfying $\text{IP}_{q}=\text{EA}_{q+1}$ of localized defects (C$_\text{B}$) in \textit{h}-BN yield gaps in good agreement with experiment \cite{expgap} and $G_0W_0$@PBE. Note that the spin-orbit coupling (which was not included in this calculation) will lower the band gap of MoS$_2$ by 0.1 eV~\cite{Alidoust2014}, which will bring our PBE0($\alpha$) and $G_0W_0$  results in even better agreement with experimental electronic band gap.}
      \label{table:gap}
\end{center}
\end{table*}


\subsection{Generalized Koopmans' Condition for the Fraction of Fock Exchange in 2D Materials}
After we obtained reliable CTLs relative to vacuum, we focused on how to calculate accurate band edge positions and band gaps of 2D materials in order to determine defect ionization energies. Using the GW approximation, we obtained an accurate quasiparticle band gap (indirect at T$\rightarrow$M) 6.01 eV for bulk h-BN (Table \ref{table:gap}), in excellent agreement with the experimental fundamental electronic gap 6.08 $\pm$ 0.015 \cite{expgap}. Meanwhile, our GW results for the band gap of bulk \textit{h}-BN (7.01 eV) and MoS$_2$ (2.82 eV), agree well with previously reported experimental values (see Table~\ref{table:gap}).
Nonetheless, GW is still computationally too demanding for defects' screening and difficult to obtain forces and optimize geometries. Therefore, the development of computationally affordable methods such as accurate non-empirical hybrid functionals for 2D materials is strongly desired. 

\begin{figure}[b]
\begin{center}
\includegraphics[keepaspectratio=true,width=0.9\linewidth]{exx_eig-eps-converted-to.pdf}
\caption{The IP at $q=0$ and the EA at $q=+1$ for the defects C$_\text{B}$ , C$_\text{N}$  and V$_\text{N}$C$_\text{B}$  in monolayer h-BN as a function of the fraction of Fock exchange $\alpha$ for PBE0($\alpha$). The predicted exchange constant ($\alpha = 0.409$, 0.41 and 0.382, respectively) is the corresponding crossing point where $\text{EA}_{q+1}=\text{IP}_{q}$.} \label{fig:ipea}
\end{center}
\end{figure}

\begin{figure}[b]
\begin{center}
\includegraphics[keepaspectratio=true,scale=0.7]{p0a_gw2-eps-converted-to.pdf}
\caption{Comparing computed band gaps of \textit{h}-BN (monolayer, bilayer, trilayer, bulk), graphane (CH), and MoS$_{2}$ with PBE0($\alpha$)\text{} versus those computed with $G_0W_0$$@$PBE . Overall we find that our PBE0($\alpha$)\text{} results agree very well with $G_0W_0$ , yielding a MAE of 0.11 eV. The grey diagonal area highlights $\pm 0.25$ eV agreement.} \label{fig:p0a}
\end{center}
\end{figure}

\begin{table}[H]
\centering
\begin{tabular}{p{2.2cm}p{1.8cm}p{1.8cm}p{1.8cm}}
\hline\hline
 System & Defect & $\alpha$ & Gap \\
\hline
ML BN 
 & C$_\text{B}$    & 0.409 & 7.344 \\
 & C$_\text{N}$    & 0.418 & 7.401 \\
 & V$_\text{N}$C$_\text{B}$  & 0.382 & 7.174 \\
BL BN   
 & C$_\text{B}$    & 0.347 & 7.075 \\
 & C$_\text{N}$    & 0.351 & 7.101 \\
 & V$_\text{N}$C$_\text{B}$  & 0.318 & 6.892 \\
TL BN  
 & C$_\text{B}$    & 0.324 & 6.994 \\
 & C$_\text{N}$    & 0.326 & 7.007 \\
 & V$_\text{N}$C$_\text{B}$  & 0.298 & 6.801 \\
Bulk BN  
 & C$_\text{B}$    & 0.225 & 6.071 \\
 & C$_\text{N}$    & 0.227 & 6.087 \\
 & V$_\text{N}$C$_\text{B}$  & 0.178 & 5.684 \\
Graphane
 & B$_\text{C}$ & 0.467 & 6.503 \\
 & N$_\text{C}$ & 0.473 & 6.541 \\
 \hline\hline
\end{tabular}
\caption{\label{table:ipea} Predicted fraction of Fock exchange for use in the PBE0($\alpha$)\text{} functional based on the IP$_{q}=$EA$_{q+1}$ condition. Note that for \textit{h}-BN the corresponding MAE compared with $G_0W_0$  results are 0.14 eV, 0.16 eV and 0.18 eV for C$_\text{B}$ , C$_\text{N}$  and V$_\text{N}$C$_\text{B}$ , respectively.}
\end{table}

The generalized Koopmans' condition has been most commonly used to determine the appropriate fraction of Fock exchange ($\alpha$) for molecules and molecular crystals \cite{Atalla2016, Thomas2014,Abrams2012,sai2011,cohen2007,Pinheiro2015,Jie2016}. 
This principle has been successfully extended to defects in bulk semiconductors \cite{miceli} through the enforcement of this condition (\textit{i.e.} $\text{EA}_{q+1}=\text{IP}_{q}$) on defects in bulk semiconductors to obtain $\alpha$ and in turn predict accurate electronic structure of the corresponding pristine bulk systems.
The fundamental assumption is that the optimized $\alpha$ depends on the long range screening of the system and not on the nature of the probe defects. This condition is also valid for deep defects in 2D materials, where defect wavefunctions are well localized like molecular orbitals in the supercells, and their  contribution to dielectric screening is negligible compared to the crystal environment.  
Another advantage of applying this condition to 2D systems is that both $\text{EA}_{q+1}$ and $\text{IP}_{q}$ can be exactly referenced to vacuum.
In order to validate the applicability of the generalized Koopmans' condition to 2D materials, we used the defect C$_\text{B}$  as a probe to determine $\alpha$ for \textit{h}-BN (B$_\text{C}$ (boron substitution of carbon) for graphane and V$_\text{S}$ (sulfur vacancy) for MoS$_2$). This method gives $\alpha$ of 
0.409, 0.347, 0.324, 0.225 for monolayer, bilayer, trilayer and bulk \textit{h}-BN, respectively as shown in Table~\ref{table:ipea}. Note that the $\alpha$ value 0.225 for bulk h-BN, agrees well with the predicted $\alpha$ from the inverse of high frequency dielectric constant ($\alpha = 1/\varepsilon_{\infty} \approx 0.2$) \cite{dielectric1}, which supports the assumption that long-range screening determines $\alpha$. We also investigated other defects C$_{\text{N}}$ and {V$_\text{N}$C$_\text{B}$} as probes of $\alpha$ as shown in Fig.~\ref{fig:ipea} (their corresponding electronic structure can be found in SI). 

Interestingly, we found that IP$_q$ and EA$_{q+1}$ from Kohn-Sham eigenvalues varied linearly with $\alpha$. 
Fig. \ref{fig:ipea} shows this linearity for three defects in monolayer \textit{h}-BN, and three defects predict similar $\alpha$, which justifies the insensitivity of $\alpha$ to the explicit defect. Note that defects with localized wavefunctions such as atomic substitutions (C$_\text{B}$ ) determine more accurate $\alpha$ compared with less localized defects such as {V$_\text{N}$C$_\text{B}$}, because the former's contribution to dielectric screening is negligible and the screening is mostly determined by the crystal. 
It is also notable that the slopes of  IP$_q$ and EA$_{q+1}$ are opposite but nearly equal, explaining why the average of IP$_q$ and EA$_{q+1}$ as $\varepsilon^{QP}$ for CTL in Eq. (\ref{eq:avg2}) works well (as shown in Table \ref{table:path}).

Most commonly, two-dimensional systems are synthesized with a few layers of the material, therefore understanding the effect of increasing thickness is essential to connect with realistic experiments. As such, we have computed the band gaps of monolayer, bilayer, trilayer and bulk \textit{h}-BN, as well as graphane and MoS$_2$, with several hybrid functionals including HSE, PBE0, B3PW and PBE0($\alpha$)\text{} (with $\alpha$ predicted earlier), and $G_0W_0$ @PBE as shown in Table \ref{table:gap} (the experimental photoemission gaps are also shown). As anticipated, PBE strongly underestimated monolayer \textit{h}-BN band gap: 4.71 eV with a direct transition at the K point. With any level of theory beyond PBE, monolayer \textit{h}-BN is predicted to have a larger, indirect gap from K to $\Gamma$. In accordance with quantum confinement, we observed that the band gaps of \textit{h}-BN obtained at B3PW,  PBE0($\alpha$) , and $G_0W_0$  show a sharp increase at ultrathin BN (monolayer to trilayer) compared to bulk BN (in agreement with a previous study \cite{wickramaratne2018monolayer}). However, HSE and PBE0 provide almost the same band gaps between ultrathin and bulk BN. This is because there is a severe change in the dielectric screening from monolayer to bulk, and a different portion of Fock exchange must be instilled. 

Using PBE0($\alpha$)\text{} we obtained results consistent with quantum confinement and in best agreement with our $G_0W_0$  calculations with a MAE of 0.11 eV (Fig.~\ref{fig:p0a}). In addition, the B3PW functional \cite{becke,bill} provided a more accurate bulk BN band gap than PBE0 and HSE but still underestimated the band gaps of ultrathin BN. 
In graphane, PBE0($\alpha$)\text{} (with a predicted exact exchange of 0.473) leads to a gap of 6.54 eV, in agreement with $G_0W_0$ @PBE, 6.41 eV. Coincidentally, MoS$_2$ has a predicted exchange of 0.250 (the default of the PBE0 scheme) and yields a gap of 2.85 eV in great agreement with $G_0W_0$ @PBE, 2.82 eV. 
Therefore, the direct/indirect transitions and magnitude of the gaps from bulk to monolayer are provided accurately solely with PBE0($\alpha$)\text{} and $G_0W_0$. 
In brief, the results shown in Table \ref{table:gap} validate our method for determining accurate fundamental band gaps for 2D materials from first-principles. We note that calculated band edge positions relative to vacuum are also similar at PBE0($\alpha$)\text{} and $G_0W_0$  as shown in Fig.~\ref{fig:ctl2} and SI.
\begin{figure*}
\begin{center}
\begin{picture}(500,180)
\put(0,10){\subfloat[\hspace{2mm}PBE\hspace{1mm}]{\includegraphics[keepaspectratio=true,scale=0.7]{cbctl_pbe-eps-converted-to.pdf}}}
\put(120,10){\subfloat[\hspace{2mm}HSE\hspace{1mm}]{\includegraphics[keepaspectratio=true,scale=0.7]{cbctl_hse-eps-converted-to.pdf}}}
\put(240,10){\subfloat[\hspace{0mm}PBE0($\alpha$)]{\includegraphics[keepaspectratio=true,scale=0.7]{cbctl_pbe0ie-eps-converted-to.pdf}}}
\put(360,10){\subfloat[\hspace{2mm}$G_0W_0$ \hspace{1mm}]{\includegraphics[keepaspectratio=true,scale=0.7]{cbctl_gw-eps-converted-to.pdf}}}
\put(305,10){\includegraphics[keepaspectratio=true,scale=0.3]{white.png}}
\put(240,10){\includegraphics[keepaspectratio=true,scale=0.7]{cbctl_pbe0ie-eps-converted-to.pdf}}
\put(185,10){\includegraphics[keepaspectratio=true,scale=0.3]{white.png}}
\put(120,10){\includegraphics[keepaspectratio=true,scale=0.7]{cbctl_hse-eps-converted-to.pdf}}
\put(65,10){\includegraphics[keepaspectratio=true,scale=0.3]{white.png}}
\put(0,10){\includegraphics[keepaspectratio=true,scale=0.7]{cbctl_pbe-eps-converted-to.pdf}}
\end{picture}
\caption{Charge transition level C$_\text{B}$  (+1/0) in \textit{h}-BN with varying number of layers at different levels of theory. The red, green, blue lines represent the band edges and defect charge transition levels for monolayer, bilayer and trilayer  \textit{h}-BN, respectively. The dashed blue line represents the charge transition level of the defect in the middle layer of the trilayer \textit{h}-BN. Defect charge transition levels gradually become shallower with lower ionization energies while increasing the number of layers (ionization energies are written adjacent to arrows from the CTL to CBM). Note that the defect CTLs are very similar relative to vacuum between different levels of theory.}  \label{fig:ctl2}
\end{center}
\end{figure*}


\subsection{Defect Ionization Energies in 2D Materials}
Finally, we can obtain the defect ionization energies based on the methods proposed earlier for CTLs and band edge positions relative to vacuum.
For example, CTLs and ionization energies for C$_\text{B}$  in \textit{h}-BN computed at PBE, HSE, PBE0($\alpha$)\text{} and $G_0W_0$ levels of theory as a function of number of layers are shown in Fig.~\ref{fig:ctl2}. Consistent with the findings in Table \ref{table:path}, CTLs changed less than 0.1 eV across different theoretical methods relative to vacuum. Interestingly, no clear trend and only small difference have been found in the band edge positions of \textit{h}-BN from monolayer to trilayer. These results illustrate that one just needs to correct the band edge positions of pristine \textit{h}-BN with PBE0($\alpha$)\text{} or $G_0W_0$ , and use CTLs determined from DFT with semi-local functionals, then the difference of the two yields accurate defect ionization energies. On another note, we found there is a clear monotonic decrease in the ionization energies of defects in \textit{h}-BN with increasing number of layers (the ionization energy is placed adjacent to the corresponding arrow in Fig.~\ref{fig:ctl2}). As shown in Fig.~\ref{fig:ie}, the ionization energy of the C$_\text{B}$  defect in \textit{h}-BN is lowered with increasing number of layers. This effect can be understood as a result of increased dielectric screening with more layers of \textit{h}-BN, and is consistent with the effect of dielectric environments on the ionization energies of MoS$_2$ \cite{mos2}. Due to accurate band gaps and the insensitivity of charge transition levels to the level of theory as discussed earlier, the PBE0($\alpha$)\text{} results agree well with $G_0W_0$ @PBE for the ionization energies of C$_\text{B}$  in \textit{h}-BN with different number of layers(blue and green points in Fig.~\ref{fig:ie}).
	
\begin{figure}[t]
\begin{center}
\includegraphics[keepaspectratio=true,scale=0.9]{ie-eps-converted-to.pdf}
\caption{Ionization energies of C$_\text{B}$  in \textit{h}-BN with varying number of layers. It is observed that ionization energies decrease monotonically with increasing number of layers. Note that PBE0($\alpha$)\text{} and $G_0W_0$  give results in excellent agreement.} \label{fig:ie}
\end{center}
\end{figure}


\section{Conclusion}

In summary, we have established fundamental principles to reliably and efficiently compute ionization energies for defects in 2D materials. 
Specifically, band edge positions of the pristine systems should be computed with our proposed PBE0($\alpha$)\text{} hybrid functional or GW approximations. We note that we obtained a fraction of Fock exchange $\alpha$ from 0.2 (bulk \textit{h}-BN) to 0.4 (monolayer \textit{h}-BN) by enforcing the generalized Koopmans' condition.
Meanwhile, we have demonstrated the insensitivity of CTL's computed by various functional choices as well as $G_0W_0$ , when the CTL's are referenced to vacuum. Therefore, the defect CTL may be obtained reliably by standard DFT with semi-local functional, when it is calculated relative to vacuum.  
We have successfully applied the proposed methods for a variety of defects from monolayer to trilayer \textit{h}-BN, graphane and MoS$_2$. The combination of these methods will allow for reliable prediction and validation of defect ionization energies in two-dimensional materials, which can be potentially used to identify the chemical composition of defects in 2D materials through comparing with experimental measurements. We also demonstrated that defect ionization energies decreased with increasing number of layers in \textit{h}-BN, due to enlarged dielectric screening. Our findings in this work suggest efficient and accurate methods to compute defect ionization energies and electronic structures in 2D materials, which can be applied to screening new promising defects for quantum information and optoelectronic applications.





\section*{Acknowledgment}
We thank Giulia Galli and Chris Van de Walle for  helpful discussions. This work is supported by NSF award DMR-1760260 and DMR-1747426, and Hellman Fellowship. T.J.S. acknowledges financial support from a GAANN fellowship. M.G. is supported by MICCoM, as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. This research used resources of the Center for Functional Nanomaterials, which is a US DOE Office of Science Facility, and the Scientific Data and Computing center, a component of the Computational Science Initiative, at Brookhaven National Laboratory under Contract No. DE-SC0012704, the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract No. DEAC02-05CH11231, the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1548562 \cite{xsede}, and the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357.

