\section{Introduction}
\label{I}
A quantum phase transition is a zero temperature transition of a quantum many body system driven by a non-commuting term of the quantum Hamiltonian
which is associated with a diverging length as well as a diverging time scale \cite{sachdev99,chakrabarti96}. In recent years, a
 plethora of studies are being carried out which attempt to bridge a connection
between quantum phase transition and quantum information theory \cite{nielsen00,vedral07}. For example, information theoritic measures like  entanglement, quantum
 fidelity \cite{zanardi06,zhou08,gu10,gritsev10}, decoherence \cite{zurek91,haroche98,zurek03,joos03}
and quantum discord \cite{sarandy09,dillenschneider08}, etc.,
are being studied close to the quantum critical point (QCP). These measures not only capture the singularities associated
with the QCP but also show distinct scaling relations which characterizes it. There have also been numerous studies on decoherence (or loss of phase information) in
 a quantum critical system which is closely connected to the LE to be discussed in this work; understanding decoherence is essential for successful achievement
of the quantum computation.

To study the LE in a quantum critical environment, we make resort to the
central spin model \cite{zanardi} in which a central spin $S$ is coupled globally to an environmental spin model $E$ (which in this case is the
two dimensional Kitaev model).
 The LE (with the $E$ in some ground state $|\psi_0\rangle$) is given by
$$ L(t) = |\langle \psi_0| e^{iH_0t} e^{-i(H_0 + \delta H_\delta)t} |\psi_0\rangle|^2. $$
Here the $H_0$ and $H_0+\delta H_\delta$ are the two Hamiltonians with which the ground state $|\psi_0\rangle$ evolves, where the term $\delta H_\delta$
 arises due to the coupling of the $E$ with the $S$.
It has been established that the LE shows a decay near
 the critical point of the $E$ with a decay rate that marks the universality associated with the QCP of $E$ \cite{zanardi, zhang09, rossini07, sharma12}. Also
the LE shows collapse
and revival as a function of time when the $E$ is at the QCP.

The proposed work is organized in the following way: Sec.\ref{I} presents the model Hamiltonian, the phase diagram and discussion about
the AQCP. In sec.\ref{II}, we describe
the general calculation of the LE and in the subsequent subsections we study the scaling of the
short time decay close to the AQCP and its collapse and revival with time.

\section{Model, Phase Diagram and Anisotropic Quantum Critical Point (AQCP)}
\label{II}
The Hamiltonian of the Kitaev model on a honeycomb lattice is given by
\begin{equation}
H=\sum_{j+l=even}\left(J_{1}\sigma^{x}_{j,l}\sigma^{x}_{j+1,l}+J_{2}\sigma^{y}_{j-1,l}\sigma^{y}_{j,l}+J_{3}\sigma^{z}_{j,l}\sigma^{z}_{j,l+1}\right)
\label{ham1}
\end{equation}
where $j$ and $l$ signify the column and row indices respectively of the honeycomb lattice  while  $J_1$, $J_2$ and $J_3$ are
coupling parameters for the three bonds (see Fig.~(\ref{Fig1})) \cite{kitaev06, sengupta08}; and $\sigma^{\alpha}_{j,l}$,  are the Pauli spin matrices
with $\alpha$($=x$, $y$ and $z$), denoting the spin component.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=7.7cm]{kit_le1.eps}
\end{center}
\caption{(Color online) Kitaev model on a honeycomb lattice with $\vec{M_1}$ and $\vec{M_2}$ being spanning
vectors of the lattice and $J_1$, $J_2$ and $J_3$, the coupling on the three bonds [\onlinecite{kitaev06}].}
\label{Fig1}
\end{figure}

We will assume the parameters $J_1$, $J_2$ and $J_3$ are all positive and
confine our analysis on the plane $J_1+J_2+J_3=4$, since only the ratio of the coupling parameters appear in the subsequent calcuations. The
most exciting property of this model is that even in two dimensions  it can be exactly solved using Jordan-Wigner (JW)
 transformation\cite{lieb61,barouch70,kogut79,bunder99,kitaev06,chen,feng}
in terms of Majorana fermions  given by
\begin{eqnarray}
&& a_{j,l}=\left(\prod_{i=-\infty}^{j-1}\sigma_{i,l}^z\right)\sigma_{j,l}^y \hspace{2mm}\text{for even $j+l$,}\nonumber\\
&& a'_{j,l}=\left(\prod_{i=-\infty}^{j-1}\sigma_{i,l}^z\right)\sigma_{j,l}^x\hspace{2mm}\text{for even $j+l$,}\nonumber\\
&& b_{j,l}=\left(\prod_{i=-\infty}^{j-1}\sigma_{i,l}^z\right)\sigma_{j,l}^x\hspace{2mm}\text{for odd $j+l$,}\nonumber\\
&& b'_{j,l}=\left(\prod_{i=-\infty}^{j-1}\sigma_{i,l}^z\right)\sigma_{j,l}^y\hspace{2mm}\text{for odd $j+l$.}
\label{majoranafermions}
\end{eqnarray}
Here,  $a_{j,l}$, $a'_{j,l}$, $b_{j,l}$ and $b'_{j,l}$ are all Majorana fermion operators, they obey the
relations $a^{\dagger}_{j,l}=a_{j,l}, b^{\dagger}_{j,l}=b_{j,l},
\{a_{j,l},a_{m,n}\}= \{b_{j,l},b_{m,n}\}= 2\delta_{j,m}\delta_{l,n}$ and $\{a_{j,l},b_{m,l}\}=0$. One can now change the lattice site indices $\left(j,l\right)$
 of honeycomb lattice to a 2-dimensional vector $\vec{n}$ , where $\vec{n}=\sqrt{3}\hat{i}n_1+
\left(\frac{\sqrt{3}}{2}\hat{i}+\frac{3}{2}\hat{j}\right)n_2$ which labels the midpoints of the vertical bonds of the honeycomb lattice. Here $n_1$ and $n_2$
take all integer values
 so that the vectors $\vec{n}$ form a triangular lattice. The Majorana fermions $a_{\vec{n}}$ and $b_{\vec{n}}$ are placed at the top and bottom
sites respectively
of the bond labeled by $\vec{n}$. The whole lattice is spanned by the vectors $\vec{M_1}= \frac{\sqrt{3}}{2}\hat{i}-\frac{3}{2}\hat{j}$
 and $\vec{M_2}= \frac{\sqrt{3}}{2}\hat{i}+\frac{3}{2}\hat{j}$, see Fig.~(\ref{Fig1}).

Under the  transformation to Majorana fermions as defined in Eqs.~(\ref{majoranafermions}), Hamiltonian (\ref{ham1}) takes the form
\begin{equation}
H=i\sum_{\vec{n}}\left(J_1b_{\vec{n}}a_{\vec{n}-\vec{M_1}}+J_2b_{\vec{n}}a_{\vec{n}+\vec{M_2}}+J_3D_{\vec{n}}b_{\vec{n}}a_{\vec{n}}\right),
\label{ham2}
\end{equation}
where $D_{\vec{n}}$=$ib'_{\vec{n}}a'_{\vec{n}}$. These $D_{\vec{n}}$ operators have eigenvalues $\pm1$ independently for each $\vec{n}$ and commute with each other
and also with $H$ which makes the Kitaev model exactly solvable. Since $D_{\vec{n}}$ is a constant of motion one can use one of the eigenvalues $\pm1$ for each
$\vec{n}$ in the Hamiltonian. The ground state of the model corresponds to $D_{\vec{n}}=1~ \forall~ \vec{n}$~ \cite{kitaev06}.
With  $D_{\vec{n}}=1$,  we can easily diagonalize the Hamiltonian (\ref{ham2}) quadratic in Majorana fermions.

The Fourier transform of the Majorana fermions can be defined as
\begin{equation}
a_{\vec{n}}= \sqrt{\frac{4}{N}}\sum_{\vec{k}}\left(a_{\vec{k}}e^{i\vec{k}.\vec{n}}+a_{\vec{k}}^{\dagger} e^{-i\vec{k}.\vec{n}}\right),
\label{ftransform}
\end{equation}
similarly $b_{\vec{n}}$ also has same Fourier transform relation. The $a_{\vec{k}}$'s and $b_{\vec{k}}$'s are Dirac fermions which follow the fermionic
anti-commutation
relations. Here, $N$ is the total number of sites and $N/2$ is the number of unit cells. In the above sum given in Eq.~(\ref{ftransform}), $\vec{k}$ is
 extended over half of the
Brillouin zone of the hexagonal lattice due to Majorana nature of the fermions \cite{sengupta08}. We recall that the full Brillouin
zone on the reciprocal lattice represents a rhombus with
vertices $\left(k_x,k_y\right)$ = $\left(\pm2\pi\sqrt{3},0\right)$ and $\left(0,\pm2\pi/3\right)$. In the momentum space the Hamiltonian (\ref{ham2}) takes the form
$H=\sum_{\vec{k}}\psi^{\dagger}_{\vec{k}}H_{\vec{k}}\psi_{\vec{k}}$ where $\psi^{\dagger}_{\vec{k}}=\left(a^{\dagger}_{\vec{k}},b^{\dagger}_{\vec{k}}\right)$ and
the reduced $2 \times 2$ Hamiltonian
 $H_{\vec{k}}$, can be expressed  in terms of Pauli matrices as

\begin{eqnarray}
&& H_{\vec{k}}= \alpha_{\vec{k}}\sigma^1 + \beta_{\vec{k}}\sigma^2,\nonumber\\
\text{where} \hspace{2mm} && \alpha_{\vec{k}} = 2[J_1\sin(\vec{k}.\vec{M_1})-J_2\sin(\vec{k}.\vec{M_2})],\nonumber\\
\text{and}\hspace{1mm} && \beta_{\vec{k}} = 2[J_3+J_1\cos(\vec{k}.\vec{M_1})+J_2\cos(\vec{k}.\vec{M_2})].
\label{ham3}
\end{eqnarray}
The eigenenergies of the $H_k$ are given by
\begin{equation}
 E^{\pm}_{\vec{k}}= \pm \sqrt{\alpha^2_{\vec{k}}+\beta^2_{\vec{k}}}.
\label{energy1}
\end{equation}
This energy spectrum corresponds to two energy bands; it is noteworthy that for $|J_1-J_2|\leq J_3\leq(J_1+J_2)$, the band
 gap $\Delta_{\vec{k}}= E^+_{\vec{k}} - E^-_{\vec{k}}$ vanishes
for some particular $\vec{k}$ modes leading to the gapless phase of the Kitaev model.
The phase diagram
of the model is shown in  an equilateral triangle  satisfying the relation $J_1+J_2+J_3=4$ and
$J_1,J_2,J_3>0$ (see Fig.~(\ref{Fig2})); one can easily show that the whole phase is divided into
three gapped phases, separated by a gapless phase (inner equilateral triangle) which is
bounded by gapless critical lines
$J_1=J_2+J_3$, $J_2=J_3+J_1$ and $J_3=J_1+J_2$.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=7.7cm]{kit_le2.eps}
\end{center}
\caption{(Color online) Figure shows phase diagram of the Kitaev model, satisfying  $J_1+J_2+J_3=4$. The inner equilateral triangle corresponds to the
gapless phase in which the coupling parameters satisfies the relations $J_1\leq J_2+J_3$, $J_2\leq J_3+J_1$
and $J_3\leq J_1+J_2$. Along the three paths I, II and III $J_3$ is varied, so as to
 study the LE. The path I, II and III are defined by the equations $J_1=J_2$, $J_1=J2+1$ and $J_1+J_3=4$ respectively.}

\label{Fig2}
\end{figure}

 On the critical line
$J_3=J_1+J_2$ energy gap goes to zero for the four $\vec{k}$ modes given by $(k_x,k_y)$= $(\pm2\pi/\sqrt{3},0)$ and $(\pm2\pi/3,0)$ which are the
four corner points of Brillouin zone. One can now expand $\alpha_{\vec{k}}$ and $\beta_{\vec{k}}$ around one of the critical modes for $J_3=J_1+J_2$, in the form
\begin{eqnarray}
&& \alpha_{\vec{k}}=  \sqrt{3}(J_2-J_1)k_x + 3(J_1+J_2)k_y,\nonumber\\
&& \beta_{\vec{k}}= \frac{3}{4}(J_1+J_2)k^2_x + \frac{9}{4}(J_1+J_2)k^2_y + \frac{3\sqrt{3}}{2}(J_2-J_1)k_xk_y,\nonumber\\
\label{alphabeta1}
\end{eqnarray}
where $k_x$ and $k_y$ are the deviations from the above mentioned critical modes.  We note that  $\alpha_{\vec{k}}$ varies linearly
and  $\beta_{\vec{k}}$ varies  quadratically
 in $k_x$ and $k_y$.  The point $J_{3c}=2J_1$ (where $J_1=J_2$) denoted by A in (Fig.~(\ref{Fig2})) needs to be checked
 carefully. This is an AQCP \cite{hikichi} with energy dispersion
 $E_{\vec{k}}\sim k^2_x$ along $k_x$
 ($k_y=0$) and $E_{\vec{k}}\sim k_y$ along $k_y$ ($k_x=0$).

The corresponding dynamical exponents are given by $z_{\perp}=1$ and $z_{\parallel}=2$, respectively.

For $J_1\neq J_2$, $J_{3c}=J_1+J_2$ is  also an AQCP that can be shown using a rotation to
a new coordinate system \begin{eqnarray}
&& k_1= \sqrt{3}(J_2-J_1)k_x + 3(J_1+J_2)k_y,\nonumber\\
&& k_2= 3(J_1+J_2)k_x - \sqrt{3}(J_2-J_1)k_y, \label{k1k2}
\end{eqnarray}
in which $\alpha_{\vec k}$ and $\beta_{\vec k}$ take the form
\begin{eqnarray}
&& \alpha_{\vec{k}}= k_1,\nonumber\\
&& \beta_{\vec{k}}= c_1k^2_1 + c_2k^2_2 + c_3k_1k_2,
\end{eqnarray}
where  $c_1 = 9(J_1+J_2)(4J^2_1+4J^2_2+J_1J_2)$,
$c_2 = 27J_1J_2(J_1+J_2)$, and
$ c_3 = 18\sqrt{3}J_1J_2(J_2-J_1)$ .
Therefore, for a general AQCP $(J_1\neq J_2)$ the dispersion will vary linearly and quadratically along $\hat{k_1}$ and $\hat{k_2}$ directions, respectively, with two
dynamical exponents $z_1=1$ and $z_2=2$.

\section{Qubit Coupled To $J_{3}$ Term Of the Kitaev Hamiltonian}
\label{III}
In this section we will provide a general calculation of the LE considering Kitaev model on a honeycomb lattice as an environment ($E$) that is coupled to a
central spin-$\frac{1}{2}$ ($S$).
 We shall denote the ground state and excited state of the central spin $S$
by $\ket{g}$ and $\ket{e}$ respectively. $S$ is coupled to $J_3$ term of $E$ Hamiltonian only when the central spin is in the excited state
$\ket{e}$. Therefore the composite Hamiltonian takes the form
\begin{eqnarray}
H_T\left(J_3,\delta\right) &=& \sum_{j+l=even}(J_{1}\sigma^{x}_{j,l}\sigma^{x}_{j+1,l}+ J_{2}\sigma^{y}_{j-1,l}\sigma^{y}_{j,l}\nonumber\\
&+& J_{3}\sigma^{z}_{j,l}\sigma^{z}_{j,l+1}+\delta\ket{e}\bra{e}\sigma^{z}_{j,l}\sigma^{z}_{j,l+1}),
\label{ham4}
\end{eqnarray}
where $\delta$ is the coupling strength of $S$ to $E$. We shall work in the limit of $\delta\rightarrow0$.

We consider  that the $S$ is initially in a generalized state $\ket{\phi(0)}_S=c_g\ket{g}+c_e\ket{e}$ (with the coefficients satisfying the condition
$|c_g|^2+|c_e|^2=1$), and the $E$ is  initially  in the ground state $\ket{\varphi(J_3,0)}$. The evolution of the environmental spin model
splits into two branches, given by $\ket{\varphi(J_3,t)}= \exp(-iH(J_3)t)\ket{\varphi(J_3,0)}$ and
$\ket{\varphi(J_3+\delta,t)}= \exp(-iH(J_3+\delta)t)\ket{\varphi(J_3,0)}$;  the evolution of $\ket{\varphi(J_3,t)}$ is driven by the
Hamiltonian $H(J_3)=H_T(J_3,0)$ (when the $S$ is in the ground state and hence there is no $\delta$ term present in the Hamiltonian),
whereas $\ket{\varphi(J_3+\delta,t)}$ evolves with $H(J_3+\delta)=H_T(J_3,0)+V_e$, where $V_e= \delta\sum_{j+l=\text{even}}\sigma^z_{j,l}\sigma^z_{j,l+1}$,
 is the effective potential arising due to the coupling between $S$ and $E$.  The wave function of the composite system
at a time $t$ is given by
\begin{equation}
\ket{\psi(t)}= c_g\ket{g}\otimes\ket{\varphi(J_3,t)} + c_e\ket{e}\otimes\ket{\varphi(J_3+\delta,t)}.
\label{twf}
\end{equation}
As a result the LE is given by
\begin{eqnarray}
 L\left(J_3,t\right) &=& |{\langle\varphi(J_3,t)|\varphi(J_3+\delta,t)\rangle}|^2,\nonumber\\
 &=& |{\langle \varphi(J_3,0)|\exp(-iH(J_3+\delta)t)|\varphi(J_3,0)\rangle}|^2.\label{le}
\end{eqnarray}
Here, we have exploited the fact that the $\ket{\varphi(J_3,0)}$ is an eigenstate of the Hamiltonian $H(J_3)$.

Following Fourier transformation and Bogoliubov transformation the diagonalized form of the Hamiltonian (\ref{ham1}) is given by
\begin{equation}
H(J_3) = \sum_{\vec{k}}[-\varepsilon_{\vec{k}}(J_3)A^{\dagger}_{\vec{k}}A_{\vec{k}} + \varepsilon_{\vec{k}}(J_3)B^{\dagger}_{\vec{k}}B_{\vec{k}}],
\label{ham6}
\end{equation}
where the $A_{\vec{k}}$'s and $B_{\vec{k}}$'s are Bogoliubov fermionic operators defined as
\begin{eqnarray}
&& A_{\vec{k}}=\frac{1}{\sqrt{2}}[a_{\vec{k}} - e^{-i\theta_{\vec{k}}}b_{\vec{k}}],
\hspace{1mm} B_{\vec{k}}=\frac{1}{\sqrt{2}}[a_{\vec{k}} + e^{-i\theta_{\vec{k}}}b_{\vec{k}}],\nonumber\\
\text{with}\hspace{1mm}&& e^{i\theta_{\vec{k}}}= \frac{\alpha_{\vec{k}}+i\beta_{\vec{k}}}{\sqrt{\alpha^2_{\vec{k}}+\beta^2_{\vec{k}}}},\label{boperators}
\end{eqnarray}
 and the
energy spectrum is given by (see Eqs.(\ref{ham3}) and (\ref{energy1}))
\begin{equation}
\varepsilon_{\vec{k}}(J_3)= \sqrt{\alpha^2_{\vec{k}} + \beta^2_{\vec{k}}}
\hspace{2mm}\text{and}\hspace{2mm} \varepsilon_{\vec{k}}(J_3+\delta)= \sqrt{\alpha^2_{\vec{k}} + \beta'^2_{\vec{k}}},
\label{energy2}
\end{equation}
where $\alpha_{\vec{k}}$ and $\beta_{\vec{k}}$ are defined in Eq.(\ref{ham3}), and $\beta'_{\vec{k}}$ corresponds to
 the value with $J_3+\delta$ instead of $J_3$.

The complete ground state of $H(J_3)$ can be written in the form (see Ref. [\onlinecite{victor12}] for details),
\begin{equation}
\ket{\varphi(J_3,0)} = \prod_{\vec{k}}[\frac{1}{2}(a^{\dagger}_{\vec{k}}-e^{i\theta_{\vec{k}}}b^{\dagger}_{\vec{k}})\hspace{2mm}(a'^{\dagger}_{\vec{k}}+ib'^{\dagger}_{\vec{k}})]\ket{\Phi}
\label{gstate}
\end{equation}
where $\vec{k}$ runs over half of the Brillouin zone of the hexagonal lattice. Following mathematical
steps identical to those in [\onlinecite{zanardi,sharma12}], it can be shown that Eq.(\ref{gstate}) leads to the expression for the LE  given by
\begin{equation}
L(J_{3},t)=\prod_{\vec{k}}L_{\vec{k}} = \prod_{\vec{k}}[1-\sin^2(2\phi_{\vec{k}})\sin^2(\varepsilon_{\vec{k}}(J_3+\delta) t)], \label{le1}
\end{equation}
where, $\tan\theta_{\vec{k}}(J_3+\delta)= \alpha_{\vec{k}}/\beta'_{\vec{k}}$ and $\phi_{\vec{k}}= [\theta_{\vec{k}}(J_3) - \theta_{\vec{k}}(J_3+\delta)]/2$.
The expression for  LE closely resembles that of  the case when the transverse Ising chain is chosen to be the environment \cite{zanardi}. For
numerical analysis of Eq.(\ref{le1}), we shall use $k_x$ and $k_y$ in terms of two
independent variables $v_1$ and
$v_2$, with $0\leq v_1,v_2\leq 1 $. The $k_x$ and $k_y$ are given by \cite{sengupta08}
\begin{equation}
k_x=\frac{2\pi}{\sqrt{3}}(v_1+v_2-1), \hspace{3mm} k_y=\frac{2\pi}{3}(v_2-v_1),
\label{kxky}
\end{equation}
which span the rhombus uniformly. Avoiding the corner points of the Brillouin zone ( where the LE results in a  zero value), we
vary $v_1$ and $v_2$ from $1/(2N)$ to $1-1/(2N)$ in steps of $1/N$, where $N$ is the system size \cite{victor12} and consider only the half of the Brillouin
 zone using the condition $(v_1, v_2)\geq 1$.

The LE is calculated numerically as a function of $J_3$ using Eq.(\ref{le1}) and it shows dip
at all critical points. To illustrate this, we choose three paths along which the interaction $J_3$ is varied.  In the first case  $J_3$ is varied along the
path $J_1=J_2$ (path `I' in Fig.~(\ref{Fig2})) so that the model  enters from the gapped phase to  the gapless phase (extending in the region  $J_3\in[0,2]$)

crossing the AQCP (point `A' in Fig.~(\ref{Fig2})) at
 $J_3=2-\delta$.  The LE shows a sharp dip at point A and there is a revival with a small magnitude which again decays at the end
 point B, $J_3=0$ (see Fig.~(\ref{Fig3})). Now surprising result shows up when the path is so chosen (path `II' in Fig.~(\ref{Fig2}), given
 by the equation $J_1=J_2+1$)
 that the system enters the gapless phase through an AQCP with $J_1 \neq J_2$, (denoted by `P' in Fig.~(\ref{Fig2})). The LE
 shows a sharp dip at the point P and stays close to its minimum value (with a small revival as
 observed in path I) throughout the gapless phase and
 again shows a rise when the system exits the gapless phase through the point Q.
 In contrary, for the case when $J_3$  is changed along the line $J_1+J_3=4$, $J_2 =0$ (path `III' in Fig.~(\ref{Fig2})), one
 observes only a single drop in the LE near $J_3=2-\delta$ (see Fig.~(\ref{Fig3}), inset (b)); this is associated with the critical point of the
 one-dimensional Kitaev model.

 In the next section
we will study the scaling of the short time behaviour of LE close to these critical points and the collapse and reviaval
of LE with time when the $E$ is right at the critical point.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=7.7cm]{kit_le3.eps}
\end{center}
\caption{ LE as a function of parameter $J_3$ ($J_3$ is varied along path 'I') shows a sharp dip at point A ($J_3=2-\delta$) and after a small revival in the gapless phase it again decays at ponit B,
 $J_3=0$ (see Fig.~(\ref{Fig2})) with $N_x=N_y=200$, $\delta=0.01$ and
$t=10$. Inset (a) shows the variation in LE when the parameter $J_3$ is varied along the path II ($J_1=J_2+1$) in the phase diagram for
$N_x=N_y=200$, $\delta=0.01$ and
$t=10$ clearly showing a sharp dip at point P ($J_3=2-\delta$) and again rise at point Q ($J_3=1-\delta$).
 Inset (b) marks the dip in LE when $J_3$ is varied along the path III ($J_1+J_3=4$), for this case $N=400$, $\delta=0.01$ and
$t=10$ so that $E$ realise the change in the behavior at $J_3=2-\delta$. Details of these three cases  is provided in the subsection (\ref{a}), (\ref{b})
 and (\ref{c}) respectively. }
\label{Fig3}
\end{figure}

\subsection{Path I: Anisotropic Quantum Critical Point ($J_1=J_2$)}
\label{a}
As discussed in Sec.\ref{II}, $J_{3c}=2J_1$ is an AQCP with critical exponents $\nu_{\perp}=z_{\perp}=1$ along $\hat j$ direction and
$\nu_{\parallel}=1/2,z_{\parallel}=2$ along
$\hat i$ direction. At this point energy gap vanishes for the three critical modes  given by $(2\pi/\sqrt{3},0)$ and $(0,\pm 2\pi/3)$ in half
of the Brillouin zone. Now
we will study the short time behavior of the LE (in Eq.~(\ref{le1})) close to the AQCP.
We define a cutoff frequency $K_c=\left(K_{x,c},K_{y,c}\right)$ such that modes up to this cut-off only are
considered to calculate the decay of LE at short time close to the AQCP. Then the LE is given by
\begin{equation}
L_c\left(J_3,t\right)=\prod_{k_x,k_y>0}^{K_c}L_{\vec{k}}.
\label{lc1}
\end{equation}
We define the quantity $\mathcal{S}(J_3,t)$, such that $\mathcal{S}\left(J_3,t\right)=\ln L_c\equiv-\sum_{k_x,k_y>0}^{K_c}|\ln L_{\vec{k}}|$. Expanding
around one of the critical mode upto the cut-off, we get
 $\sin^2\varepsilon_k(J_3+\delta)t \approx 4\left(J_3+\delta-2J_1\right)^2t^2$ and $\sin^2\left(2\phi_k\right)\approx9J_1^2k_y^2\delta^2/\left(J_3-2J_1\right)^2
\left(J_3+\delta-2J_1\right)^2$ therefore we obtain,
\begin{equation}
\mathcal{S}\left(J_3,t\right) \approx -\frac{36\mathcal{E}\left(K_c\right)J_1^2\delta^2t^2}{\left(J_3-2J_1\right)^2},
\label{sj3}
\end{equation}
where $\mathcal{E}\left(K_c\right)$=$4\pi^2N_c\left(N_c+1\right)\left(2N_c+1\right)/54N_y^2$ and $N_c$ is a integer nearest to $3N_yK_c/2\pi$.
We therefore find an exponential decay of the LE in the early time limit given by
\begin{equation}
L_c\left(J_3,t\right)\approx \exp\left(-\gamma t^2\right)
\label{lc2}
\end{equation}
where $\gamma=36\mathcal{E}\left(K_c\right)J_1^2\delta^2/\left(J_3-2J_1\right)^2$.
The anisotropic nature of the quantum critical point is reflected in the fact that $\gamma$ scales as $1/N_y^2$ and is independent of $N_x$.
Further using the expression of $L_c\left(J_3,t\right)$, one can easily observe that it is invariant under the transformation $N_y\rightarrow N_y \alpha, \delta\rightarrow \delta/\alpha$ and $t\rightarrow t\alpha$, where $\alpha$ is some integer.

Now we fix $J_1=J_2=1$ and $J_3=2-\delta$ (point `A' in Fig.~(\ref{Fig2})) and observe collapse and revival of LE with time (presented in the Fig.~(\ref{Fig4})). The time period of collapse and revival is proportional to $N_y$,
and is unaffected by the changes in $N_x$; this confirms the scaling result of the decay rate $\gamma$ for the short time limit near the AQCP discussed above.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8.7cm, height=6.5cm]{kit_le4.eps}
\end{center}
\caption{ (Color online) The collapse and revival of LE with t at the AQCP (point `A' in Fig.~(\ref{Fig2})) for $J_1=J_2=1$, $\delta=0.01$ and
$J_3=2-\delta$, keeping $N_x (=100)$, fixed and varying $N_y$ verifies
 the scaling relations satisfied by $N_y$, $\delta$ and $t$ as discussed in the text. The inset shows that the quasiperiod of the collapse and revival is independent of $N_x$. }
\label{Fig4}
\end{figure}

\subsection{Path II: Anisotropic Quantum Critical Point $\left(J_1 \neq J_2\right)$}
\label{b}
It has been shown that the point P in Fig.~(\ref{Fig2}) ($J_1\neq J_2$, $J_{3,c}=J_1+J_2 $)
is an AQCP which can be seen by choosing directions
$\hat{k_1}= \sqrt{3}\left(J_2-J_1\right)\hat{i} + 3\left(J_1+J_2\right)\hat{j} $ (see Sec \ref{II}) and $\hat{k_2}$, perpendicular
 to $\hat{k_1}$ \cite{victor12}.
The critical exponents associated with this critical point are given by $\nu_1=z_1=1$, and $\nu_2=1/2,z_2=2$ along $\hat k_1$ and $\hat k_2$ directions, respectively.
To calculate the early time scaling in a similar spirit as in the previous section,
we expand Eq.~(\ref{le1}) near one of the critical modes up to the cut-off $K_c$ to  obtain
$\sin^2\left(\varepsilon_k(J_3+\delta)\right)t \approx 4\left(J_3+\delta-J_1-J_2\right)^2t^2$ and $\sin^2\left(2\phi_k\right)\approx k_1^2\delta^2/
4\left(J_3-J_1-J_2\right)^2\left(J_3+\delta-J_1-J_2\right)^2$. In the short time limit, the LE becomes
\begin{equation}
L_c\left(J_3,t\right) \approx \exp\left(-\gamma t^2\right)
\label{lc3}
\end{equation}
where, $\gamma=\delta^2\mathcal{E}\left(K_c\right)/\left(J_3-J_1-J_2\right)^2$, $\mathcal{E}\left(K_c\right)= 8\pi^2J_2^2N_c\left(N_c+1\right)\left(2N_c+1\right)/3N^2$ and
$N_c$ is an integer nearest to $NK_c/4\pi J_2$.

In fact comparing with the previous section \ref{a}, one can see that in this case $k_1$ (instead of $k_y$) appears in
the expression  of the LE in the short-time limit.
 Further, from
 Eq. (\ref{lc3}) and the expression of $\gamma$ one observes that $L_c\left(J_3,t\right)$ is invariant under the
transformation $N_x=N_y=N\rightarrow N\alpha,\delta\rightarrow\delta/\alpha$ and
$t\rightarrow t\alpha$, with $\alpha$ being some integer
which is also observed in the collapse and revival behavior (see Fig.~(\ref{Fig5})).
\begin{figure}[ht]
\begin{center}
\includegraphics[width=7.7cm]{kit_le5.eps}
\end{center}
\caption{ (Color online) Variation of LE with $t$ at the AQCP P ($J_1=3/2,J_2=1/2$ and $J_3=2-\delta$) shows collapse
and revival with diffrent $N_x=N_y=N$ and $\delta=0.001$. }
\label{Fig5}
\end{figure}

\subsection{Path III: One-dimensional Quantum Critical Point $\left(J_2=0\right)$}
\label{c}
As mentioned already, along the line $J_1+J_3=4$,  ($J_2=0$), the two dimensional spin model reduces to an equivalent one dimensional  spin chain with energy gap vanishing at $J_1=J_3$ for
$k_c=\pi$ and the corresponding dynamical exponent being $z=1$. We shall now expand $\sin\left(\varepsilon_k(J_3,\delta)t\right)$ and $\sin\left(2\phi_k\right)$ around the critical mode $k_c$ to analyze the short time
decay of LE, resulting into $\sin^2\left(\varepsilon_k(J_3+\delta)t\right)\approx 4\left(J_3+\delta-J_1\right)^2t^2$ and
$\sin^2\left(2\phi_k\right)\approx J_1^2k^2\delta^2/\left(J_3-J_1\right)^2\left(J_3+\delta-J_1\right)^2$. The LE  hence takes the from
\begin{equation}
L_c\left(J_3,t\right)\approx \exp\left(-\gamma t^2\right)
\label{lc4}
\end{equation}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=7.7cm]{kit_le6.eps}
\end{center}
\caption{ (Color online) The LE as function of time at the QCP (point `R' in Fig.~(\ref{Fig2})) $J_1=2$ and $J_3=2-\delta$ with different N and $\delta=0.01$, verifying the analytical
scaling with $N$, $\delta$ and $t$. }
\label{Fig6}
\end{figure}
where, $\gamma=4J_1^2\delta^2\mathcal{E}\left(K_c\right)/\left(J_3-J_1\right)^2$ and $\mathcal{E}\left(K_c\right)= 4\pi^2N_c\left(N_c+1\right)\left(2N_c+1\right)/6N^2$ , $N_c$ is an integer nearest to
$NK_c/2\pi$. From Eq.~(\ref{lc4}), we find  that the LE shows a similar scaling relation as is
expected for a one-dimensional chain with $z=1$ \cite{zanardi};  this is also confirmed by studying the
collapse and revival of LE (see Fig.~(\ref{Fig6})).

\section{Conclusion}
In this paper we study a variant of the central spin model in which a central spin (qubit) is globally coupled to an environment which is chosen to be a
two-dimensional  Kitaev model on a honeycomb lattice through the interaction term $J_3$. Using the exact solvability of the Kitaev model,
we have derived  an exact expression of the LE when the interaction $J_3$ is varied in a way such that
the system enters the gapless phase crossing the AQCP of the phase diagram. However, the
behavior of the LE as a function of $J_3$ depends upon the path along which $J_3$ is varied. In the
case when the AQCP, Q (with $J_1 \neq J_2$ see Fig.~(\ref{Fig2})) is crossed, one observes a complete revival of the echo when
the system exits the gapless phase to re-enter the gapped phase; this is in contrast
to the case $J_1 =J_2$. For the case of $J_2=0$ there is only one sharp dip at the critical point $J_3=2-\delta$ which is associated
with the QCP of the one-dimensional Kitaev model.

The early time scaling behavior for both the paths I and II close to the AQCP bear the signature
of the fact that the gapless phase is entered crossing an AQCP with different exponents along
different spatial directions. This is also confirmed by studying the collapse and the revival of
the LE as a function of time. However, one does not observe a perfect collapse and revival (except
for the equivalent one dimensional case); this may be because of the proximity to a gapless phase. The quasi-period of collapse and revival
 in all cases scale with the system size as $N^z$.  The case with $J_2=0$ reflects the fact that the system is essentially
one-dimensional in this limit.  It is straightforward to relate these results to the decoherence
of the central spin close to a critical point.

{This study of LE can  be verified experimentally as presented by Zhang $et~al$ \cite{zhang09}; they measure the LE as an indicator
 of quantum criticality for a
one-dimensional quantum Ising model with an antiferromagnetic interaction using NMR quantum
simulators. In this experiment, they prepare the ground state of the Hamiltonian (using the gate sequences) which need not
be the true ground state but could be a state that
approximates the ground state of the system well, and then measure the LE for finite number of spins. Similar experiments can be realized with
the approximate ground state of the Kitaev model. Also, since the Kitaev model can also be realized using an optical lattice\cite{duan03, sen08}
(where the couplings can be separately tuned with
the help of different microwave radiations), there exists a possibility of verifying these results in an optical lattice also.}

{It should be noted that  in a recent work, Pollmann $et~al$ \cite{pollmann10}, have studied the problem of the LE in a transverse
Ising spin chain in the presence
of a longitudinal field; more precisely they calculated
the magnitude of the overlap between the final state reached following a slow quench across the
QCP and its time evolved counterpart at time $t$ (generated following the time evolution with the final Hamiltonian). They
observe a cusp-like minimum in the echo as a function of time in the
limit when the spin chain is integrable. However, this behavior is smeared in the non-integrable
case (with non-zero longitudinal field) thus providing a probe for integrable versus non-integrable
behavior. In the present paper, we  however deal with an equilibrium situation in which the spin
chain is not quenched across the QCP, and observe the collapse and revival only at the QCP.}

\begin{center}
\bf{Acknowledgements}
\end{center}
We acknowledge  Amit Dutta, Victor Mukherjee and Aavishkar Patel for helpful discussions and comments. AR acknowledges B. K.
Chakrabarti for valuable discussions and thanks IIT Kanpur for financial support during this work. SS thanks CSIR, New Delhi for Junior Research Fellowship.


\section{Introduction}

Despite the positive attributes of organic materials for next-generation photovoltaics, the effective deployment of organic photovoltaic (OPV) cells poses fundamental problems at the molecular scale. Besides enhancing the durability of organic compounds, the central requisite to the mass-market viability of OPV modules is to raise their power conversion efficiency (PCE) beyond 10-15\%. \cite{HachmannOlivares-Amaya2011} Comparatively, the efficiency of current OPV architectures (Fig.~\ref{OSC}) barely exceeds 9-10\% PCE under standard test conditions. \cite{DouYou2012}  Nevertheless, the margin for improvement is vast due to the exceptional chemical versatility of organic materials, \cite{ScharberMuhlbacher2006} and accurate approaches to predict the charge-transfer and optical properties of new compounds could allow for extensive screening of promising semiconducting organics. \cite{SokolovAtahan-Evrenk2011}

\begin{figure}[ht!]
\includegraphics[width=8.5cm]{OSC.eps}
\caption{Schematics of an organic solar cell showing a bilayered donor-acceptor heterojunction (top). Bilayered, nanopatterned, and mixed-phase heterojunctions with typical power conversion efficiencies (bottom).
\label{OSC}}
\end{figure}

However, conventional density-functional theory approximations in their static (DFT) and time-dependent (TDDFT) formulations \cite{HohenbergKohn1964, KohnSham1965, RungeGross1984} suffer from severe limitations in predicting the electronic structure of semiconducting materials, \cite{CohenMori-Sanchez2008} precluding the quantitative and even qualitative elucidation of relevant photovoltaic mechanisms. Notably, conventional DFT and TDDFT approximations tend to systematically destabilize occupied states and overstabilize unoccupied states in organic and organometallic complexes, leading in particular to the incorrect description of the electronic properties (charged excitations) and optical properties (neutral excitations) of donor-acceptor dyads, extended polymer molecules, and heavy-metal sensitizing dyes. \cite{DreuwWeisman2003, FaassenBoeij2002, BlaseAttaccalite2011a, HimmetogluMarchenko2012} To overcome DFT limitations, one privileged route has been to resort to many-body approaches, namely, to many-body perturbation theory approximations such as GW. \cite{Hedin1965,OnidaReining2002} The GW method, while accurate in predicting electronic spectra, is much more expensive than DFT although considerable progress has been achieved in reducing the cost of GW calculations. \cite{UmariStenuit2009,UmariStenuit2010, BlaseAttaccalite2011,BlaseAttaccalite2011a,FaberAttaccalite2011,QianUmari2011}

In parallel to many-body perturbation theory, less expensive orbital-dependent density-functional theories (OD-DFTs) represent promising alternatives. \cite{KummelKronik2008} At present, the most widely used OD-DFT methods are hybrid density-functional theory  approximations. \cite{Becke1993} Hybrid functionals in their linear-admixture or range-separated forms have been shown to improve upon conventional DFT in predicting electronic properties. \cite{KummelKronik2008, VarsanoAndrea-Marini2008, BroqvistAlkauskas2009, SteinEisenberg2010} Alternatively, self-interaction corrections to density-functional theory  approximations \cite{PerdewZunger1981} represent a second category of OD-DFT methods that are now rapidly growing in recognition due to their improved accuracy and lower computational cost. In self-interaction-corrected approaches, the total energy of the system is expressed explicitly in terms of individual orbital densities to rectify nonphysical errors inherent in orbital-independent DFT. Successful applications of advanced self-interaction corrections have appeared with an accuracy potentially approaching that of many-body methods in predicting electronic levels and electrical responses. \cite{BaerNeuhauser2005, KorzdorferMundt2008, DaboFerretti2010, FerrettiDabo2012}

In this work, we highlight the predictive accuracy of the recently developed self-interaction correction based upon Koopmans' theorem (Koopmans-compliant OD-DFT) \cite{DaboFerretti2010, FerrettiDabo2012} in capturing the electronic levels of OPV compounds. Electronic levels represent reliable PCE indicators \cite{ScharberMuhlbacher2006} and have successfully served as fundamental inputs for first-principles combinatorial screening of organic semiconductors. \cite{HachmannOlivares-Amaya2011, SokolovAtahan-Evrenk2011} We demonstrate that Koopmans-compliant OD-DFT is apt at describing donor levels within 0.1-0.4 eV and acceptor levels within 0.2-0.6 eV relative to experiment, which is comparable to the precision of the GW method. In its simplest formulation, the method can be trivially implemented in conventional electronic-structure codes with an algorithmic cost potentially lower than that of hybrid-DFT functionals as it does not require evaluating exchange terms. Furthermore, the method holds promise for accurate time-dependent extensions based upon the improved description of the underlying electronic spectrum.

This work is organized as follows. First, we underscore the significance of Koopmans' theorem in understanding the charge-transfer behavior of donor-acceptor complexes. Second, we outline existing electronic-structure approaches that aim at enforcing Koopmans' theorem. Third, we present the Koopmans-compliant OD-DFT method to correct the nonphysical tendency of conventional DFT approximations to destabilize occupied states and overstabilize unoccupied states. Last, we demonstrate the efficiency and accuracy of the method for families of organic compounds of interest to photovoltaics.

\section{Method}

\label{Methodology}

\subsection{Koopmans' theorem}

\begin{figure}[ht!]
\includegraphics[width=7.5cm]{InfiniteSeparation.eps}
\caption{Energy of a positively charged pentacene-fullerene dyad in the infinite separation limit (black) and energies of the isolated fractionally charged pentacene P$^q$ (light blue) and fullerene C$^{1-q}_{60}$ molecules (dark blue) as a function of the transferred charge $q$ within the LSDA and Koopmans-compliant descriptions. (Energies are given relative to the fully ionized individual molecules.)
\label{InfiniteSeparation}}
\end{figure}

Koopmans-compliant functionals aim at imposing Koopmans' theorem in DFT approximations. Koopmans' theorem enables one to equate orbital energies (that is, the expectation values of the effective Hamiltonian) with total energy differences, which correspond to withdrawing an electron from a stationary electron state.

 The relevance of Koopmans' theorem on the accuracy of effective-potential theories in predicting electronic spectra has been recognized in several theoretical and computational studies. \cite{HachmannOlivares-Amaya2011, SteinEisenberg2010, BaerLivshits2010, FerreiraMarques2011, AlkauskasPasquarello2011, JarlborgBarbiellini2011, CohenMori-Sanchez2012, KronikStein2012} In addition, Koopmans' theorem is the central condition to correctly describe electron transfer within donor-acceptor pairs. To illustrate this fact, we consider a positively charged pentacene-fullerene dyad in the limit of infinite intermolecular separation where electronic interactions between the two molecules can be neglected.

Figure \ref{InfiniteSeparation} depicts the energies of the isolated pentacene oligomer and fullerene cluster. The total energy of the infinitely separated molecular pair as a function of the transferred charge $q$ is also reported. Within the local spin-density approximation (LSDA) [Fig.~\ref{InfiniteSeparation}(a)], the three curves exhibit a marked nonlinear dependence on $q$, reflecting the fact that LSDA violates Koopmans' theorem. In other words, the derivative of the energy $\textrm{d}E(\textrm{X}^q)/\textrm{d}q$ of an individual fractionally charged cation X$^q$, which corresponds to the  opposite energy of its highest occupied orbital, departs significantly from the finite energy difference $E(\textrm{X}^+)-E(\textrm{X})$, which corresponds to the ionization potential of the neutral molecule X. More precisely, the dependence of the energies is found to be convex, causing the dyad to be most stable at a fractional charge $q_\textrm{min}$ close to $\frac 23$. This observation is in qualitative contradiction to the expectation that the charge should fully localize onto the strongest electron acceptor (namely, fullerene). In contrast to this nonphysical behavior, Koopmans-compliant energy curves are linear [Fig.~\ref{InfiniteSeparation}(b)] as a result of Koopmans' condition $\textrm{d}E(\textrm{X}^q)/\textrm{d}q=E(\textrm{X}^+)-E(\textrm{X})$, leading to the expected stabilization of the fully transferred charge ($q_\textrm{min}=1$).

Similar observations would be made for other donor-acceptor complexes. Moreover, at finite separation where electronic coupling between the two molecules must be taken into account, LSDA would also significantly underestimate the transferred charge. These results highlight the central importance of imposing Koopmans' theorem on the individual subsystems to correctly describe charge transfer in donor-acceptor pairs.

\subsection{Existing methods}

\label{ExistingMethods}

At present, there exist different DFT-based methods that aim at imposing Koopmans' theorem. \cite{Dabo2008, SalznerBaer2009, LanyZunger2009, SteinEisenberg2010, BaerLivshits2010, LanyZunger2010, Refaely-AbramsonBaer2011} The self-interaction correction of Perdew and Zunger \cite{PerdewZunger1981} can be regarded as one of the earliest approaches to restore the physical interpretation of orbital energies as ionization energies. In this approach, the Kohn-Sham density functional $E^{\rm KS}$ is augmented with corrections that remove the self-interaction of the electronic wave functions in the one-electron limit. In explicit terms, the Perdew-Zunger (PZ) correction to the Kohn-Sham DFT functional reads $E^{\rm PZ}=E^{\rm KS}-\sum_{i\sigma} E^{\rm KS}_{\rm Hxc}[\rho_{i\sigma}]$, where $\rho_{i\sigma}$ stands for the density of the one-electron state $\psi_{i\sigma}$ and $E^{\rm KS}_{\rm Hxc}$ denotes the Hartree plus exchange-correlation component of the Kohn-Sham functional. The PZ approach is more expensive than local and semilocal DFT approximations but appreciably less costly than hybrid-DFT approximations. It rectifies the tendency of local and semilocal DFT to destabilize occupied orbitals in atoms but overcorrects the error for molecules, leading to an overestimation of ionization energies of up to 1-2 eV. Also, the PZ correction leaves the energies of empty levels unchanged, meaning that acceptor levels remain largely underestimated (that is, excessively negative). Furthermore, the lack of balance in the PZ correction of occupied and unoccupied orbitals leads to the underestimation of the total energies of molecules and the underestimation of interatomic distances.

A solution to the lack of balance in the PZ correction has been proposed by Baer et al. in the form of a range-separated hybrid functional [the Baer-Neuhauser-Livshits (BNL) functional] in which long-range interactions are treated at the explicit-exchange level, whereas short-range interactions are described at the semilocal density-functional level. \cite{BaerNeuhauser2005,Livshitsbaer2007} In order to impose Koopmans' theorem within BNL, the original empirical parameter $\gamma$ that determines the range of the electrostatic separation ($1/r={\rm erf}(\gamma r)/r + {\rm erfc}(\gamma r)/r$) is optimally tuned to impose the agreement between orbital energies and ionization energies for the frontier orbitals of the system, \cite{SalznerBaer2009,BaerLivshits2010} leading to considerable improvement in predicting electronic and optical spectra. \cite{SteinEisenberg2010, Refaely-AbramsonBaer2011} In terms of computational cost, BNL calculations are as expensive as hybrid-DFT calculations.

A recently proposed approach that targets minimal computational cost is the method of Lany and Zunger, \cite{LanyZunger2009, LanyZunger2010} which consists of rectifying deviations from Koopmans' theorem through on-site projection terms that shift the electronic levels in order to impose the correspondence between orbital energies and ionization energies. Because it only relies on the calculation of on-site occupations, the Lany-Zunger method does not cause any noticeable increase in computational cost relative to the underlying DFT approximation. However, it requires to preselect atomic orbitals to construct the on-site projectors, making its application less systematic than the PZ correction and optimally tuned BNL method.

\subsection{Present method}

The Koopmans-compliant method that we have introduced in Ref.~\onlinecite{DaboFerretti2010} provides a promising alternative to existing approaches. The method is systematic and allies conceptual simplicity, computational performance, and predictive accuracy, thereby permitting the precise and efficient electronic-structure description of a wide variety of compounds.

The motivation for introducing the Koopmans-compliant method originates from the practical observation that determining ionization potentials (IPs) of molecules as differences of ground-state energies $I_N=E_{N-1}-E_N$ (where $I_N$ and $E_N$ denote the IP and total energy of the system of $N$ electrons) yields predictions in very good agreement with experiment within conventional DFT approximations, namely, the LSDA approach and semilocal generalized-gradient approximations (GGAs) --- in the literature, the procedure that consists in evaluating ionization energies from differences of DFT total energies is commonly known as the $\Delta$SCF method. \cite{KowalczykYost2011} For electron affinities (EAs), $\Delta$SCF predictions $A_N=E_N-E_{N+1}$ are also found to be in agreement with experiment with the central caveat that $\Delta$SCF predictions for EAs are only possible insofar as the system can accommodate the added electron within the approximation used; this condition on the stability of anionic states is necessary to calculate the ground-state energy of negatively charged systems. In practical terms, failure to stabilize anions is frequent within approximate DFT \cite{KimSim2011} and arises from spurious self-interaction whereby an electron can interact with itself through the nonphysical contribution from its own charge density to the effective Kohn-Sham Hamiltonian. \cite{PerdewZunger1981}

Therefore, in order to obtain reliable predictions for orbital levels (especially, for frontier donor and acceptor levels that determine the charge-transfer properties of the system), it is crucial to satisfy the following conditions: (a) the correspondence between the effective orbital levels and $\Delta$SCF total-energy differences must be imposed; (b) the established accuracy of $\Delta$SCF energies must be preserved; and (c) the calculations must circumvent the inability of approximate functionals to stabilize negatively charged systems. Fulfilling these requirements is the object of the non-Koopmans OD-DFT correction that is described below.

\begin{figure}[ht!]
\includegraphics[width=8.5cm]{CarbonAdiabatic.eps}
\caption{Diagram of the non-Koopmans energy $\Pi^{\rm u}_{\rm 2p}$, which quantifies  departure from linearity for the adiabatic ionization curve of the 2p state of carbon within LSDA.
\label{CarbonAdiabatic}}
\end{figure}

In presenting the non-Koopmans correction, it is important to distinguish between electronically adiabatic (i.e., unrelaxed) ionization processes, whereby orbitals are kept frozen upon ionizing the system, and diabatic (i.e., relaxed) ionization processes, whereby orbitals are allowed to rearrange self-consistently. For adiabatic ionization, the correspondence between orbital energies and total energies is referred to as the restricted Koopmans' theorem, whereas for diabatic ionization, it is termed the generalized Koopmans' theorem.\cite{SzaboOstlund1989} In what follows, we consider the adiabatic case before addressing the diabatic case.

Bearing in mind this distinction, the starting point of the method consists in defining an adiabatic measure of the deviation from Koopmans' theorem in terms of an orbital-dependent energy functional. We thus consider for the moment an unrelaxed process that corresponds to depleting the stationary state $\psi_{i\sigma}$. For this process, the energy measure that monitors deviations from the restricted Koopmans' theorem can be written as the non-Koopmans energy
\begin{equation}
\Pi^{\rm u}_{i\sigma}(f)=\int_0^{f_{i\sigma}} {\rm d}f' [\epsilon^{\rm u}_{i\sigma}(f) - \epsilon^{\rm u}_{i\sigma}(f')],
\label{NKEa}
\end{equation}
first introduced by Perdew and Zunger in discussing the nonlinear behavior of the DFT total energy upon ionization. \cite{PerdewZunger1981} In Eq.~\eqref{NKEa}, the superscript notation  $\{\cdot\}^{\rm u}$ indicates that the orbitals are kept unrelaxed along the ionization path, $f_{i\sigma}$ denotes the initial occupation of the ionized orbital, and $\epsilon^{\rm u}_{i\sigma}$ stands for the expectation value of the effective Kohn-Sham Hamiltonian for the ionized state $\psi_{i\sigma}$ along the unrelaxed curve. In other words, $\epsilon^{\rm u}_{i\sigma}$ equals the derivative of the unrelaxed total energy $E^{\rm u}_{i\sigma}$ as a function of the occupation of the ionized state (the restricted Janak theorem). Alternatively, the unrelaxed non-Koopmans energy can be written in explicit terms as
\begin{eqnarray}
\Pi^{\rm u}_{i\sigma}(f)&=&f_{i\sigma}(2f-f_{i\sigma})E_{\rm H}[n_{i\sigma}]
+E_{\rm xc}[\rho-\rho_{i\sigma}] - E_{\rm xc}[\rho] \nonumber \\
&+& \int {\rm d}{\bf r} v_{{\rm xc},\sigma}\textstyle\left({\bf r};\left[\rho+\left(f - f_{i\sigma}\right)n_{i\sigma}\right]\right)\rho_{i\sigma}({\bf r}),
\label{NKEb}
\end{eqnarray}
where $\rho_{i\sigma}({\bf r})$ stands for the orbital density, $n_{i\sigma}({\bf r})=\rho_{i\sigma}({\bf r})/\int {\rm d}{\bf r}' \rho_{i\sigma}({\bf r}')$ denotes the normalized density, $E_{\rm xc}$ and $E_{\rm H}$ represent the exchange-correlation and Hartree energy functionals, and $v_{{\rm xc},\sigma}$ is the exchange-correlation potential.

The pictorial interpretation of the non-Koopmans energy is simple. Graphically, the unrelaxed non-Koopmans energy $\Pi^{\rm u}_{i\sigma}(f)$ corresponds to the error made in approximating the ionization curve connecting the initial state $(I)$ to the final state $(F)$ by a straight line whose slope is the derivative of the total energy at the transition point $(T_f)$ (Fig.~\ref{CarbonAdiabatic}). In fact, if the adiabatic Koopmans' theorem were satisfied, that is, if the energy $\epsilon^{\rm u}_{i\sigma}(f)$ were constant and equal to the $\Delta$SCF energy difference $E^{\rm u}_{i\sigma}(f_{i \sigma})-E^{\rm u}_{i\sigma}(0)$, the ionization curve would be exactly linear and the non-Koopmans deviation $\Pi^{\rm u}_{i\sigma}(f)$ would vanish for any occupation $f$ between 0 and $f_{i\sigma}$.

With the non-Koopmans measure in hand, it becomes possible to impose Koopmans' theorem, as shown in Ref.~\onlinecite{DaboFerretti2010}. Here, for the purpose of outlining the method, the non-Koopmans-corrected (NKC) functional can be described simply as a sum of weighted orbital constraints that are added to the Kohn-Sham total energy functional to penalize deviations from Koopmans' theorem:
\begin{equation}
E^{\rm NKC} = E^{\rm KS}+{\sum_{i\sigma}} \alpha^{\rm NK}_{i\sigma} \Pi^{\rm u, KS}_{i\sigma}(f^{\rm NK}_{i\sigma}),
\label{NKC}
\end{equation}
In Eq.~\eqref{NKC}, the terms $\Pi^{\rm u, KS}_{i\sigma}$ are the adiabatic non-Koopmans deviations that correspond to the chosen Kohn-Sham approximation, the coefficients $\alpha_{i\sigma}^{\rm NK}$ represent the weights of the non-Koopmans penalties, which must be calculated self-consistently to equate orbital energies with $\Delta$SCF energy differences [condition (a)] and the coefficients $f^{\rm NK}_{i\sigma}$ stand for the reference occupations, which must be fixed to preserve the accuracy of $\Delta$SCF energies [condition (b)], as explained further below.

At this point, it is important to note that Kohn-Sham density-functional theory is a ground-state theory, which is originally not intended to predict excited-state energies. As a consequence, only the non-Koopmans energies of the highest occupied and lowest unoccupied molecular orbitals can be rigorously defined within DFT since they correspond to ionization processes that take the $N$-electron ground state into the ground states with $N-1$ and $N+1$ electrons, respectively. In contrast, the ionization of other electron orbitals ends up into non-Aufbau states whose non-Koopmans energies are in principle not defined within DFT.

As a consequence, the OD-DFT penalty sum that appears in Eq.~\eqref{NKC} should be rigorously restricted to the frontier highest occupied and lowest unoccupied states. Despite this recognized limitation, it is realized in practical computations that $\Delta$SCF based upon local and semilocal DFT functionals predict accurate excitation energies provided that orbitals do not delocalize nonphysically. \cite{KowalczykYost2011,HimmetogluMarchenko2012} These computational observations indicate that generalizing non-Koopmans penalties to the full orbital spectrum in Eq.~\eqref{NKC} should yield a Koopmans-compliant OD-DFT method able to predict electronic structures with good accuracy. In practice, it is found that this generalization provides spectral predictions in very good agreement with experimental data. \cite{DaboFerretti2010, FerrettiDabo2012}

\begin{figure}[ht!]
\flushleft (a) unrelaxed \\
\center \includegraphics[width=5cm]{CarbonAlphaAdiabatic.eps} \\
\flushleft (b) relaxed \\
\center \includegraphics[width=5cm]{CarbonAlphaDiabatic.eps}
\caption{Occupation dependence of the non-Koopmans-corrected total energy upon ionization of the 2p state of carbon as a function of the penalty coefficient $\alpha_{\rm 2p}^{\rm NK}$ in the (a) unrelaxed and (b) relaxed cases. The underlying Kohn-Sham approximation is LSDA. The reference occupation is of $f_{\rm 2p}^{\rm NK} = \frac 12$.
\label{CarbonAlpha}}
\end{figure}

To complete the presentation of the Koopmans-compliant functional, it remains to explain the determination of the penalty coefficients $\alpha_{i\sigma}^{\rm NK}$ and reference occupations $f^{\rm NK}_{i\sigma}$. The analytical determination of $\alpha_{i\sigma}^{\rm NK}$ has already been presented. \cite{DaboFerretti2010} In this analysis, it has been shown that, starting from convex-energy local and semilocal DFT functionals, there always exist a value of $\alpha_{i\sigma}^{\rm NK}$ lying between 0 and 1 that equalizes the slope of the ionization curve (i.e., the orbital energies) at the initial point $(I)$ and that at the final point $(F)$. Furthermore, the value of $\alpha_{i\sigma}^{\rm NK}$ that fulfills this condition can be calculated straightforwardly through the secant-method recursion:
\begin{equation}
\alpha_{n+1} = \alpha_n + \frac{(1-\alpha_n)\Delta^{\rm NK}_{i\sigma}(\alpha_n)}{\Delta^{\rm  NK}_{i\sigma}(\alpha_n)-\Delta^{\rm  NK}_{i\sigma}(1) },
\label{ALP}
\end{equation}
where $\Delta^{\rm NK}_{i\sigma}(\alpha)$ denotes the difference of orbital energies from $(I)$ to $(F)$ for the penalty coefficient $\alpha^{\rm NK}_{i\sigma}=\alpha$. At that precise value, the energy curve is closely linear and Koopmans' theorem is satisfied.

To illustrate this fact, we depict the influence of $\alpha_{i\sigma}^{\rm NK}$ on the linearity of the energy curve upon ionization of the 2p state of carbon in Fig.~\ref{CarbonAlpha}. Considering first the adiabatic case [Fig.~\ref{CarbonAlpha}(a)], it is seen that the linearity of the ionization energy, that is, the fulfilment of the adiabatic Koopmans' theorem, is achieved for $\alpha_{\rm 2p}^{\rm NK}=1$. This result remains valid for any adiabatic process, as demonstrated by analyzing non-Koopmans residual errors. \cite{DaboFerretti2010}

Instead, in the diabatic case, the value of $\alpha_{i\sigma}^{\rm NK}$ that satisfies the generalized Koopmans' theorem is found to be lower than 1. As a matter of fact, for relaxed ionization of the 2p state of carbon, it is seen in Fig.~\ref{CarbonAlpha}(b) that linearity is achieved for $\alpha_{\rm 2p}^{\rm NK}$ close to $\frac 45$. Precisely, the converged value of the penalty coefficient that is determined after 2 recursion steps [Eq.~\eqref{ALP}] is of 0.85. More generally, it has been shown that  $\alpha_{i\sigma}^{\rm NK}$ can be regarded as a relaxation coefficient that measures the stability of the final ionized state. Thus, for systems that consist of a lone electron around a closed shell (e.g., alkali-metal atoms), the outer-shell penalty coefficient equals 1. In contrast,  $\alpha_{i\sigma}^{\rm NK}$ becomes lower than 1 upon withdrawing an electron from a filled or partly filled shell (cf.~Fig.~9 in Ref.~\onlinecite{DaboFerretti2010}).

\begin{figure}[ht!]
\includegraphics[width=7cm]{CarbonReferenceOccupation.eps}
\caption{Occupation dependence of the non-Koopmans-corrected total energy upon relaxed ionization of the 2p state of carbon as a function of the reference occupation $f_{\rm 2p}^{\rm NK}$. The underlying Kohn-Sham approximation is LSDA. The penalty coefficient is set to be $\alpha_{\rm 2p}^{\rm NK}=0.85$, as calculated recursively from Eq.~\eqref{ALP}.
\label{CarbonReferenceOccupation}}
\end{figure}

After explaining the calculation of the penalty coefficients, we now turn our attention to the determination of the reference occupation. Based upon Slater's theorem and the sum rule satisfied by the exchange-correlation hole, it can be shown that $f_{i\sigma}^{\rm NK}$ should be set to $\frac 12$ in order to preserve the accuracy of the underlying Kohn-Sham functional in predicting total-energy differences. \cite{DaboFerretti2010} This fact is illustrated in Fig.~\ref{CarbonReferenceOccupation} where it can be seen that the Koopmans-compliant energy difference matches the Kohn-Sham $\Delta$SCF energy difference at half reference occupation for diabatic ionization of the 2p state of carbon. Therefore, at the reference occupation $f^{\rm NK}_{\rm 2p}=0.5$ and with the relaxation parameter $\alpha_{\rm 2p}^{\rm NK}=0.85$ that has been determined from Eq.~\eqref{ALP}, the Koopmans-compliant functional yields an orbital energy $\epsilon_{\rm 2p}^{\rm NK}=11.83$ eV on a par with the accuracy of the $\Delta$SCF prediction (11.72 eV) and experimental IP (11.26 eV).

These results demonstrate the possibility of rectifying the nonlinear behavior of the ionization curve and restoring Koopmans' theorem [condition (a)] while preserving the accuracy of total-energy predictions [condition (b)] within Koopmans-compliant OD-DFT. In addition, Koopmans compliance presents the central advantage of avoiding to treat negatively charge molecular states [condition (c)] for both the evaluation of non-Koopmans penalties and the calculation of orbital levels, thereby allowing us to determine accurate acceptor levels straightforwardly, as shown in Sec.~\ref{BenchmarkPredictionsAcceptor}.

Further details on the implementation of one Koopmans-compliant method (the $\alpha$NKC$_0$ computational approach) are presented in the Appendix.

\section{Results}

\label{BenchmarkPredictions}

We now assess the predictive ability of Koopmans-compliant functionals in describing the donor-acceptor properties of a range of organic molecules relevant to OPV molecular junctions, first focusing on donor levels in Sec.~\ref{BenchmarkPredictionsDonor} and then examining the problem of acceptor levels in Sec.~\ref{BenchmarkPredictionsAcceptor}.

\subsection{Donor levels}

\label{BenchmarkPredictionsDonor}

\begin{figure}[ht!]
\includegraphics[width=8.5cm]{OrganicMolecules.eps}
\caption{Representative sample of organic molecules (Ref.~\onlinecite{BlaseAttaccalite2011}) for benchmarking the performance of electronic-structure methods.
\label{OrganicMolecules}}
\end{figure}

\begin{figure}[ht!]
\includegraphics[width=8.5cm]{OrganicMoleculesDonor.eps}
\caption{LSDA and $\alpha$NKC$_0$ donor levels compared with experiment (cf.~references in Table \ref{OrganicMoleculesDonorComparison}) for the organic molecules listed in Fig.~\ref{OrganicMolecules}.
\label{OrganicMoleculesDonor}}
\end{figure}

\begin{table}
\caption{Absolute donor levels of benchmark organic molecules within LSDA, PZ, $\alpha$NKC$_0$, and GW compared with experimental vertical ionization energies (cf.~Ref.~\onlinecite{NIST2012} and references below).}
\label{OrganicMoleculesDonorComparison}
\begin{tabular*}{1.0\columnwidth}{@{\extracolsep{\fill}}lccccc}
\hline \hline \\
& \multicolumn{5}{c}{Donor levels (eV)} \\
\cline{2-6} \\
& LSDA & PZ & $\alpha$NKC$_0$ & GW\footnotemark[1] & Expt.\\
\hline
\\
benzothiadiazole & 6.40 &10.29& 9.03 & 8.56 & 9.15\footnotemark[2] \\
benzothiazole & 6.20 & 9.66 & 8.64 & 8.48 & 8.85\footnotemark[3] \\
fluorene & 5.64 & 8.69 & 7.51 & 7.64 & 7.91\footnotemark[4]\\
PTCDA & 6.32 & 9.66 & 8.12 & 7.68 & 8.2\footnotemark[5] \\
phthalocyanine & 5.21 & 8.31 & 6.54 & 6.10 & 6.41\footnotemark[6] \\
porphine & 5.24 & 8.34 & 6.66 & 6.70 & 6.9\footnotemark[7] \\
tetraphenylporphyrin & 4.89 & 8.10 & 6.21 & 6.20 & 6.39\footnotemark[8] \\
thiadiazole & 7.13 &11.73&10.10& 9.89 & 10.11\footnotemark[9]\\
thiophene & 6.02 &10.97& 8.76 & 8.63 & 8.85\footnotemark[10]\\
\\
MAD & 2.20  & 1.43 & 0.16 & 0.32 \\
RMS & 0.64 & 0.45 & 0.11 & 0.15 \\
\hline \hline
\end{tabular*}
\flushleft
\footnotemark[1]{Reference \onlinecite{BlaseAttaccalite2011}.} \\
\footnotemark[2]{Reference \onlinecite{JohnstoneMellon1973}.} \\
\footnotemark[3]{Reference \onlinecite{SalmonaFaure1975}.} \\
\footnotemark[4]{Reference \onlinecite{RuscicKovac1978}.} \\
\footnotemark[5]{Reference \onlinecite{DoriMenon2006}.} \\
\footnotemark[6]{Reference \onlinecite{Berkowitz1979}.} \\
\footnotemark[7]{Reference \onlinecite{DupuisRoberge1980}.} \\
\footnotemark[8]{Reference \onlinecite{KhandelwalRoebber1975}.} \\
\footnotemark[9]{Reference \onlinecite{PasinszkiKrebsz2010}.} \\
\footnotemark[10]{Reference \onlinecite{BajicHumski1985}.} \\
\end{table}

\begin{figure}[ht!]
\includegraphics[width=8.5cm]{AcenesDonor.eps}
\caption{LSDA and $\alpha$NKC$_0$ donor levels compared with experiment for benzene (Ref.~\onlinecite{PiancastelliKelly1987}) and oligoacenes (Ref. ~\onlinecite{BiermannSchmidt1980}).
\label{AcenesDonor}}
\end{figure}

\begin{figure}[ht!]
\includegraphics[width=8.5cm]{FullerenesDonor.eps}
\caption{LSDA and $\alpha$NKC$_0$ donor levels compared with experiment for fullerenes. Experimental vertical ionization potential are from Ref.~\onlinecite{ZimmermanEyler1991} for C$_{50}$, C$_{60}$, and C$_{80}$ (D$_{\rm 5d}$) and from Ref.~\onlinecite{LichtenbergerRempe1992} for C$_{70}$.
\label{FullerenesDonor}}
\end{figure}

The initial benchmark that is depicted in Fig.~\ref{OrganicMolecules} consists of a representative set of molecules, which has been studied in Ref.~\onlinecite{BlaseAttaccalite2011}. These molecules constitute elemental components for relevant photovoltaic organics and allow us to assess the accuracy of the self-interaction method in capturing the electronic structure of organic materials across a representative sample of chemical compositions and molecular sizes.

Koopmans-compliant calculations are carried out at the $\alpha$NKC$_0$ level. We employ conventional norm-conserving pseudopotentials to represent atomic cores with a cutoff energy of 60 Ry in expanding wave functions. All of our calculations take into account spin polarization. Employing auxiliary-function correction,\cite{LiDabo2011} a vacuum separation of 9 \AA\ is sufficient to achieve convergence of the electronic structure within 0.15 eV, even for the large PTCDA and phthalocyanine molecules. In all cases, less than two recursions of Eq.~\eqref{ALP} are needed in converging the penalty coefficient $\alpha$ of the highest occupied state of the molecule, and the same penalty weight is imposed on the full spectrum (as justified by the sensitivity analysis presented in the Appendix). Damped electronic dynamics is used in optimizing electronic degrees of freedom with the inner-loop procedure described in Ref.~\onlinecite{ParkFerretti2012} and with a convergence threshold of 10$^{-5}$ Ha on the electron-dynamics kinetic energy. Regarding the experimental literature, it is important to mention that all of the measurements that are reported here correspond to molecules isolated in the gas phase. Therefore, environmental effects arising from a surrounding solvent or from an embedding bulk need not be included in our calculations.

LSDA and $\alpha$NKC$_0$ predictions for the frontier donor levels (i.e., the absolute energy of the highest occupied orbital) of benchmark molecules are reported and compared with experiment in Fig.~\ref{OrganicMoleculesDonor}. (Here and throughout the comparative assessment, we focus on LSDA results rather than on GGA predictions due to the fact that $\alpha$NKC$_0$ is applied to LSDA; this permits to directly resolve the effective influence of $\alpha$NKC$_0$ on the underlying functional. Note that GGA orbital energies would be very close to the LSDA results on the scale of the self-interaction error.)

First, we observe a marked underestimation of absolute donor levels within LSDA. The underestimation is above 1 eV for porphyrin complexes and reaches 2.5-3 eV for strong aromatic donors. This incorrect trend is rectified by $\alpha$NKC$_0$ that reduces the error down to a few tenths of an eV. Furthermore, the accuracy of $\alpha$NKC$_0$ compares favorably with that of the original self-interaction correction (PZ)\cite{PerdewZunger1981} and GW approximation. Indeed, data reported in Table~\ref{OrganicMoleculesDonorComparison} reveal that PZ provides some improvement over LSDA, reducing the mean absolute error from 2.20 eV to 1.43 eV. Much better agreement with experiment is achieved within self-consistent GW\footnote{In these GW calculations, wave functions are kept frozen and self-consistency is performed on eigenenergies.} with a mean absolute error as low as 0.32 eV and a standard deviation of the error of 0.15 eV. Despite the very good performance of GW calculations,\footnote{It is beyond the scope of this work to enter into the delicate discussion of the accuracy of GW predictions. Notwithstanding recognized difficulties in converging GW calculations, the present comparisons provide a sufficiently precise assessment of the performance of first-principles methods to draw conclusions on their relative performance.}  the Koopmans-compliant $\alpha$NKC$_0$ method is found here to be more precise in predicting the highest donor levels with an error of 0.16 eV and a low standard deviation of 0.11 eV.

After highlighting the remarkable precision of $\alpha$NKC$_0$ across a variety of chemical compositions, we now assess its ability to capture the influence of molecular size on the electronic structure. For this complementary benchmark, we consider two important families of organic molecules, namely, acene oligomers and fullerene clusters.  We employ the same calculation parameters as those above --- with the exception of the plane-wave cutoff that can be reduced to 40 Ry without altering numerical convergence.

For oligoacenes (Fig.~\ref{AcenesDonor}), we observe that LSDA is in error of 2.8 eV for the donor level of benzene. The error decreases gradually with the length of the chain but is still as large as 2 eV for hexacene. In terms of relative errors, the LSDA underestimation fluctuates in the narrow range 29-31\%. Applying $\alpha$NKC$_0$ restores the agreement between electronic-structure calculations and experiment, reducing the error to less than 0.2 eV regardless of the length of the chain. Similarly, for fullerenes (Fig.~\ref{FullerenesDonor}), $\alpha$NKC$_0$ brings about some considerable improvement in capturing the delicate influence of cluster size. Indeed, despite the fact that the LSDA error varies in a much wider range than for acene oligomers (22-28\%) and that the dependence of the highest occupied level as a function of the molecular size is clearly nonmonotonic, the $\alpha$NK$_0$ error never exceeds 0.2 eV compared to available experimental data, providing further confirmation of the predictive precision of the Koopmans-compliant method in describing subtle electronic-structure trends.

In addition, an important advantage of the $\alpha$NKC$_0$ method lies in its moderate computational cost. In fact, converging the full electronic spectrum of an extended PTCDA molecule with a large plane-wave cutoff of 60 Ry on a conventional 16-processor machine requires half of a CPU day, whereas the same calculation within HF and hybrid-DFT is 3 times more costly using comparable {\sc cp} implementations. Indeed, as explained previously, although conventional self-interaction corrections and hybrid-DFT both scale quartically [$O(N^4)$] with the size of the electronic system, the algorithmic prefactor of self-interaction methods is lower than that of hybrid-DFT since the former do not require to evaluate pairwise exchange contributions.

\subsection{Acceptor levels}

\label{BenchmarkPredictionsAcceptor}

\begin{table}
\caption{Absolute acceptor levels of benchmark organic molecules within LSDA, PZ, $\alpha$NKC$_0$, and GW compared with available experimental data.}
\label{OrganicMoleculesAcceptorComparison}
\begin{tabular*}{1.0\columnwidth}{@{\extracolsep{\fill}}lccccc}
\hline \hline \\
&  \multicolumn{5}{c}{Acceptor levels (eV)} \\
\cline{2-6} \\
& LSDA & PZ & $\alpha$NKC$_0$ & GW\footnotemark[1] & Expt. \\
\hline
\\
benzothiadiazole & 3.52 & 3.42 & 1.07  & 0.42 & ---  \\
benzothiazole & 2.35 & 2.20 & $<$0 & $<$0 & ---   \\
fluorene & 2.05 & 2.00 & $<$0 & $<$0 & ---  \\
PTCDA & 4.80 & 4.81 & 3.19  & 2.68 & ---  \\
phthalocyanine & 3.79 & 3.82 & 2.34  & 2.07 & ---  \\
porphine & 3.28 & 3.19 & 1.53  & 1.39 & ---  \\
tetraphenylporphyrin & 3.07 & 3.04 & 1.68  & 1.49  & 1.69(10)\footnotemark[2] \\
thiadiazole & 2.95 & 2.78 & $<$0 & $<$0 & --- \\
thiophene & 1.59 & 1.38 & 0.04  & $<$0 & --- \\
\\
\hline \hline
\end{tabular*}
\flushleft
\footnotemark[1]{Reference \onlinecite{BlaseAttaccalite2011}.} \\
\footnotemark[2]{Reference \onlinecite{ChenPan1991}.}
\end{table}

After validating the $\alpha$NKC$_0$ method for occupied donor levels, we now concentrate on the prediction of unoccupied acceptor levels. As already mentioned, acceptor levels represent a central difficulty for $\Delta$SCF predictions due to the inability of conventional DFT methods to properly stabilize negatively charged states. In what follows, we demonstrate that $\alpha$NKC$_0$ provides a solution to this important methodological limitation.

In calculating acceptor levels, we use the same computational parameters as those of donor-level calculations. In particular, the $\alpha$NK$_0$ penalty weights remain unchanged. First-principles predictions for benchmark organic compounds are reported in Table \ref{OrganicMoleculesAcceptorComparison}. Our results clearly confirm the poor performance of LSDA and PZ in predicting acceptor levels; both methods overstabilize acceptor states by several eVs relative to GW due to self-attraction. In contrast, $\alpha$NKC$_0$ and GW predictions are found to be in close agreement. In quantitative terms, the discrepancy between $\alpha$NKC$_0$ and GW is as low as 0.15 eV for porphine and is at most of 0.6-0.7 eV for benzothiadiazole. It should also be noted that our comparison reveals that $\alpha$NKC$_0$ acceptor levels are always more stable than their GW counterparts. Additionally, although experimental data for benchmark molecules are scarce, comparison with the measured vertical electron affinity of tetraphenylporphyrin suggests that the predictive accuracy of $\alpha$NKC$_0$ is comparable with the accuracy of GW predictions.

\begin{table}
\caption{Absolute acceptor levels of acene chains of increasing length within LSDA, $\alpha$NKC$_0$,and GW compared with experimental data.}
\label{AcenesAcceptor}
\begin{tabular*}{1.0\columnwidth}{@{\extracolsep{\fill}}lcccc}
\hline \hline \\
& \multicolumn{4}{c}{Acceptor levels (eV)} \\
\cline{2-5}\\
& LSDA & $\alpha$NKC$_0$ & GW\footnotemark[1] & Expt. \\
\hline
\\
benzene       & 1.38 & 0.03 & --- & --- \\
naphthalene & 2.27 & 0.49 & --- & --- \\
anthracene   & 2.85 & 1.16 & 0.29 & 0.53\footnotemark[2] \\
tetracene      & 3.22 & 1.60 & 0.93 & 1.07\footnotemark[3] \\
pentacene    & 3.48 & 1.98 & 1.36 & 1.39\footnotemark[3] \\
hexacene     & 3.68 & 2.18 & ---  & --- \\
octacene      & 4.27 & 2.48 & ---  & --- \\
\\
MAD & 2.19 & 0.59 & 0.13 \\
RMS & 0.09 & 0.04 & 0.09  \\
\hline \hline
\end{tabular*}
\flushleft
\footnotemark[1]{Reference \onlinecite{BlaseAttaccalite2011}.} \\
\footnotemark[2]{Reference \onlinecite{SchiedtWeinkauf1997}.} \\
\footnotemark[3]{Reference \onlinecite{CrockerWang1993}.}
\end{table}

\begin{table}
\caption{Absolute acceptor levels of fullerene clusters of increasing size within LSDA, $\alpha$NKC$_0$, QMC (quantum Monte-Carlo), GW$_0$ (nonself-consistent), scGW$_{\rm f}$ (self-consistent with vertex corrections) compared with experimental data.}
\label{FullerenesAcceptor}
\begin{tabular*}{1.0\columnwidth}{@{\extracolsep{\fill}}lcccccc}
\hline \hline \\
& \multicolumn{6}{c}{Acceptor levels (eV)} \\
\cline{2-7}\\
& LSDA & $\alpha$NKC$_0$ & QMC\footnotemark[1] & GW$_0$\footnotemark[1] & scGW$_{\rm f}$\footnotemark[1] & Expt. \\
\hline
\\
C$_{20}$ & 4.24 & 2.19 & 1.76(11) & 3.55 & 2.36 & 2.25\footnotemark[2] \\
C$_{24}$ & 5.06 & 2.88 & 2.57(11) & 4.19 & 2.88 & --- \\
C$_{50}$ & 5.20 & 3.49 & 3.52(14) & 4.75 & 3.73 & 3.10\footnotemark[3] \\
C$_{60}$ & 4.27 & 2.64 & 2.23(19) & 3.87 & 2.98 & 2.69\footnotemark[4]  \\
C$_{70}$ & 4.10 & 2.61 & 2.46(11) & 3.98 & 2.83 & 2.76\footnotemark[5] \\
C$_{80}$ (D$_{\rm 5d}$) & 5.00 & 3.91 & 3.25(10) & 4.62 & 3.88 & 3.70\footnotemark[3] \\
C$_{80}$ (I$_{\rm h}$)    & 4.56 & 2.99 & 3.90(11) & 5.17 & 4.38 & --- \\
\\
MAD & 1.66 & 0.18 & 0.42 & 1.25 & 0.26 \\
RMS & 0.32 & 0.12 & 0.06 & 0.26 & 0.21 \\
\hline \hline
\end{tabular*}
\flushleft
\footnotemark[1]{Reference \onlinecite{TiagoKent2008}.} \\
\footnotemark[2]{Reference \onlinecite{PrinzbachWahl2006}.} \\
\footnotemark[3]{Reference \onlinecite{YangPettiette1987}.} \\
\footnotemark[4]{Reference \onlinecite{WangDing1999}.} \\
\footnotemark[5]{Reference \onlinecite{WangWoo2006}.}
\end{table}

To gain further insight into the performance of first-principles methods in predicting acceptor levels, we concentrate on  the lowest unoccupied state for acenes (Table \ref{AcenesAcceptor}) and for fullerenes (Table \ref{FullerenesAcceptor}). For weak acene acceptors, the first notable feature is the propensity of LSDA to overestimate absolute acceptor levels. This expected nonphysical tendency is rectified to a large extent by $\alpha$NKC$_0$, which reduces the mean absolute error from 2.19 to 0.59 eV relative to vertical ionization experiments. Nevertheless, it is seen that GW performs significantly better than $\alpha$NKC$_0$, bringing the error down to 0.13 eV. In the same vein, we note that the naphthalene acceptor level is predicted to be stable within $\alpha$NKC$_0$ in contradiction to experiment. In fact, among all of our benchmark calculations, acenes represent the only instance where GW performs better than $\alpha$NKC$_0$, suggesting that GW should be preferred over $\alpha$NKC$_0$ in capturing shallow acceptor levels.

Despite this fact, the scenario is entirely different for strong acceptors. Indeed, as shown in Table \ref{FullerenesAcceptor}, $\alpha$NKC$_0$ yields acceptor levels in close agreement with experimental data for fullerene clusters with a mean absolute deviation of 0.18 eV, which is significantly lower than that found for many-body quantum Monte-Carlo (0.42 eV), nonself-consistent many-body perturbation theory (1.25 eV), and self-consistent many-body perturbation theory (0.26 eV). These results suggest us to employ Koopmans-compliant $\alpha$NKC$_0$ to describe the electronic properties of the strong acceptor compounds relevant to photovoltaic applications.

Overall, the systematic comparison of electronic-structure methods offers a clear validation of $\alpha$NKC$_0$ as an accurate and efficient method to describe the donor and acceptor levels of relevant OPV compounds. This remarkable performance allows us to access the electronic structure of novel organics and provides a reliable starting point to address the charge-transfer and optical properties of donor-acceptor photovoltaic junctions.

\section{Conclusion}

In order to overcome the predictive deficiency of conventional DFT methods in capturing the donor and acceptor levels of molecular complexes, we have presented a correction procedure that enables one to eliminate the nonphysical curvature of the total energy upon withdrawing (injecting) electrons from (into) the molecular system. When applied to LSDA, the procedure yields an orbital-dependent functional that fulfills the generalized Koopmans' theorem, thereby restoring or imposing the agreement between frontier orbital levels and $\Delta$SCF energy differences, and overcoming central limitations in applying $\Delta$SCF techniques to predict acceptor levels. Specifically, the $\alpha$NKC$_0$ formulation of the Koopmans-compliant method has been applied to OPV compounds spanning a wide range of chemical compositions and molecular sizes, thereby highlighting the conceptual simplicity and computational performance of the method with an accuracy comparable to that of many-body perturbation theory.

Although the study has essentially focused on determining frontier donor and acceptor levels, the non-Koopmans corrective procedure can be straightforwardly extended to compute the full electronic spectrum of molecular systems in a single first-principles calculation.\cite{FerrettiDabo2012} For future studies, we plan to apply the Koopmans-compliant method to other relevant families of photovoltaic compounds, including metallic porphyrins and phthalocyanines, and to generalize the approach to describe strongly interacting dyads and bulk organic semiconducting materials at both the electronic and excitonic levels.

\acknowledgments

The authors are indebted to Arash Mostofi, Xavier Blase, and Kieron Burke for valuable discussions and practical suggestions. ID acknowledges partial support from the French National Research Agency through Grant ANR 12-BS04-0001 PANELS (Photovoltaics from Ab-initio Novel Electronic-structure Simulations). AF acknowledges partial support from Italian MIUR through Grant FIRB-RBFR08FOAL\_001. MC acknowledges partial support from NSF EAR 0810272 and from the NSF CAREER award DMR 1151738.


\section{Introduction}
\section{Introduction}

In the first four decades since their discovery, pulsars were almost
the exclusive domain of radio (and to a lesser extent X-ray)
astronomy. Indeed, of the $\sim$2000 known pulsars\footnote{{\tt
    http://www.atnf.csiro.au/research/pulsar/psrcat}}, the majority were
discovered in radio. Since June 2008, however, with the launch of the
LAT, on the {\it Fermi} satellite, $\gamma$-rays have become a viable means of discovering (and studying)
pulsars. More important than the number of LAT-detected pulsars
(a small fraction of the overall population), the LAT sample is subject to
different biases than the radio sample. LAT pulsars are typically nearby ($\sim$few kpc) and energetic ($\dot{E}>10^{33}$ erg s$^{-1}$),
and a large fraction are {\it radio-quiet}. Furthermore,
$\gamma$-rays, unlike the radio beams, carry a significant fraction of the rotational energy of
pulsars, thus providing a powerful probe into these extreme objects. Finally, LAT pulsars provide a crucial input into our
understanding of the overall neutron star population of the Galaxy.

The LAT is a pair conversion telescope consisting of a tracker, a
calorimeter, and a segmented anti-coincidence
detector, along with a programmable trigger and a complex data
acquisition system. By incorporating the latest advances from particle
physics\footnote{Bill Atwood was awarded the 2012 Panofsky Prize in
  Experimental Particle Physics ``For his leading work on the design,
  construction, and use of the Large Area Telescope [...]''.}, including the use of silicon-strip detectors, the LAT has
achieved a giant leap in capabilities compared to its predecessor,
EGRET. The LAT extends to higher energies ($>$300 GeV vs $\sim$10
GeV), has a larger effective area and field of view, improved angular
and energy resolutions, and much lower deadtime (\cite[Atwood et al. 2009]{Atwood09}). In its first 4
years of operations, the LAT has collected $>$200 million ``source'' class events, compared to $\sim$1.5 million
photons collected by EGRET, in its 9-year lifetime.

Regarding pulsars, the progress made by the LAT has been
spectacular. The $\gamma$-ray pulsar population has grown from 7 firm
detections (plus a few candidates) at the time of the {\it Fermi}
launch (\cite[Thompson 2008]{Thompson08}), to 117 detections in three years of LAT
survey observations (The Second {\it Fermi}-LAT Catalog of
$\gamma$-ray Pulsars, {\it in preparation}). Beyond the jump in the number of pulsars, the increased
statistics enable detailed (e.g. phase-resolved) studies on individual
pulsars, previously out of reach. The LAT has
also opened up the unexplored 10--100 GeV window. EGRET
detected a mere handful of $>$10 GeV photons from the
brightest pulsars (\cite[Thompson 2005]{Thompson05}), whereas the LAT detects significant
emission from over two dozen pulsars in this energy range and even $>$25 GeV
events from the brightest ones (e.g.  Crab, Vela, Geminga, see
\cite[Saz Parkinson 2012]{Gamma2012}).

Of the 117 LAT detections, 61 pulsars were known prior to  {\it
  Fermi}. Pulsations were obtained by folding the $\gamma$-rays with a
radio (or X-ray) timing model. Of these 61, 20 are MSPs
(e.g. \cite[Abdo et al. 2009c]{MSP}), a hitherto unknown class of strong
$\gamma$-ray emitters. More surprisingly, a large number of MSPs ($>$40)
were discovered in radio searches of LAT
unassociated sources (\cite[Ray et al. 2012]{PSC}). Most (20 so far, e.g. \cite[Ransom et al. 2011]{Ransom11}) will
likely exhibit $\gamma$-ray pulsations\footnote{Given that they are
  LAT $\gamma$-ray sources selected precisely for their {\it pulsar-like} qualities.}, once enough data and/or a precise timing model is obtained (usually from
radio observations). Finally, 36 of the 117 $\gamma$-ray
pulsars were discovered directly in blind searches of LAT data
(\cite[Abdo et al. 2009b, Saz Parkinson et al. 2010, Saz Parkinson
2011, Pletsch et al. 2012a, Pletsch et al. 2012b]{BSPI,BSPII,Sardinia,Pletsch12,Pletsch12b}). The rest of this
paper discusses the challenges, results, and prospects of these searches.

\section{Blind searches for $\gamma$-ray pulsars (compared to radio)}

Two factors make $\gamma$-ray searches for pulsars particularly
challenging, compared to radio searches. The first involves the
scarcity of events. The LAT detects $\sim$1 $\gamma$-ray per day
from a typical (bright) pulsar. This means that LAT searches for pulsars span
months to years (see Figure 1). The second complication involves the
broad, energy-dependent PSF of the instrument (from $\sim5^\circ$ at 100 MeV to $\sim0.2^\circ$ at 100 GeV,
68\% containment, normal incidence). This results in significant
source confusion, especially in high background regions like the Galactic plane, where
the diffuse $\gamma$-ray emission makes it hard to resolve individual
sources. Essentially, it is impossible to select, with certainty,
events coming from a source, so maximum likelihood techniques must be
employed. Using these techniques, it is possible to determine the
probability that each event is coming from the source of interest and
improve the sensitivity of the search by assigning a weight to each event
equal to this probability \cite[(Kerr 2011)]{Kerr11}.

\begin{figure}[htbp]
\includegraphics[height=.25\textheight]{J1745.eps}
\includegraphics[height=.25\textheight]{EMS1331_0.7_400_3_PRESTO_091005.eps}
\caption{ PRESTO (\cite[Ransom 2001]{Ransom01}) plots for two pulsars: one using radio data (left) and
  one using $\gamma$-ray data (right). {\bf Left} PSR J1745+1017, an MSP discovered in radio searches of {\it Fermi} unassociated sources
  with the Effelsberg telescope (Barr et al. 2012, MNRAS, {\it
    submitted}. Figure Credit: E. Barr). {\bf Right} PSR J1957+5033, a {\it radio-quiet} $\gamma$-ray pulsar discovered
  in a LAT blind search. Note that, although the pulse profiles and
  significances are similar, the $\gamma$-ray observation is $\sim$10,000 times
  longer, and the number of $\gamma$-ray events (582) is approximately
  one per day ($\sim$200,000 pulsar rotations!). On the other hand, unlike in radio, there is no need to
  scan in DM in a $\gamma$-ray search.}
\label{fig6}
\end{figure}

Searches over long (sparse) data sets make standard FFT techniques
impractical. To tackle this problem, \cite[Atwood et al. (2006)]{Atwood06} developed the ``time-differencing
technique''  and applied it successfully to EGRET data
(\cite[Ziegler et al. 2008]{Ziegler08}). The core of the technique involves computing the
FFT of the differences between the times of events (up to a maximum,
sliding, time window), rather than using the original time series. Using a
time window of $\sim$weeks significantly alleviates the coherence requirements of the search.

\subsection{Young pulsars}

Within weeks of the launch of {\it Fermi}, the first {\it radio-quiet}
$\gamma$-ray pulsar was discovered in the supernova remnant CTA 1(\cite[Abdo et al. 2008]{CTA1}). This was
followed by many more discoveries (e.g. \cite[Abdo et al. 2009b, Saz Parkinson et
al. 2010]{BSPI,BSPII}). Many of these early
$\gamma$-ray pulsars discovered by the LAT are coincident with old
EGRET unidentified sources. Many, like the pulsar in CTA 1, are also
coincident with known supernova remnants or pulsar wind nebulae, and
were thus long suspected of hosting pulsars (e.g. PSR J1836+5925, J2021+4026). Radio follow-up
observations of the LAT-discovered pulsars showed that most of them are {\it radio-quiet} (or extremely
radio-faint), suggesting that the $\gamma$-ray emission originates far
from the neutron star surface, resulting in broad beams. It is now
clear that these radio-quiet $\gamma$-ray pulsars represent a significant fraction
of the neutron star population of the Galaxy (see Guillemot et al. (2012), in these proceedings, for a detailed
review of the radio observations of LAT $\gamma$-ray pulsars). All 36 pulsars found in LAT blind searches to date are young
($\tau<1\times10^7$ yr) energetic ($10^{33}$ erg s$^{-1} <\dot{E}< 10^{37}$ erg s$^{-1}$) pulsars with frequencies below
$\sim$20 Hz. These pulsars often exhibit timing irregularities, in
the form of {\it timing noise} and glitches. Since most are not detected in radio, the LAT is the only instrument
capable of timing them (\cite[Ray et al. 2011]{Ray11}). While it is
possible to time these noisy (or ``glitchy'') pulsars, a good timing
model often requires many frequency derivatives and other whitening
terms, making their {\it discovery} over long data spans extremely challenging. Searches for these pulsars may only be possible
over LAT observing periods lasting $\sim$months, rather than years (see Figure 2).

\begin{figure}[b]
\begin{center}
\includegraphics[width=4.5in]{J1023_Signif_Prob.eps}
\caption{Blind-search significance of PSR J1023-5745, as a function of the cumulative
  observing time. Initially, the significance increases with time,
  but after a period of 7-10 months, the significance peaks and then
  decreases with the addition of more data.
}
\label{fig6}
\end{center}
\end{figure}

\subsection{Millisecond pulsars}

As mentioned above, the discovery of $\gamma$-ray MSPs (\cite[Abdo et
al. 2009c]{MSP}) was largely unexpected. Equally
unexpected was the large number of MSPs discovered by radio
telescopes searching LAT unassociated sources (\cite[Ray et al. 2012]{PSC}). This fact, combined
with the large fraction of {\it radio-quiet} young pulsars, raises the
question of whether large large numbers of {\it radio-quiet}
$\gamma$-ray MSPs remain to be discovered in blind searches of
LAT data. The answer, from examining the identified $\gamma$-ray
sources in the {\it almost} flux-limited {\it Fermi}-LAT Bright Source List (\cite[Abdo et al. 2009a]{BSL}),
appears to be no. Whereas at least two thirds of young $\gamma$-ray pulsars are {\it
  radio quiet}, at most one third of $\gamma$-ray MSPs are (\cite[Romani 2012]{Romani12}).

Blind searches for MSPs are vastly more complicated than searches for young pulsars. Firstly, a
majority ($\sim$80\%) of MSPs are in binary systems, and a full blind
search over unknown orbital parameters is out of reach given current
computer capabilities. Searches for isolated MSPs are possible, but
the high frequencies increase the memory and CPU
requirements. Positional uncertainties become more relevant with
increasing frequencies and the typical tolerance of a blind search
using several years of data is of order a fraction of an arc
second. This means that one must either perform a fine scan over the
LAT positions, or else identify a precise position of a plausible
counterpart, using multi-wavelength (e.g. X-ray) observations.
A number of LAT sources have been identified as promising pulsar
candidates, by virtue of their variability and spectral properties
(\cite[Ackermann et al. 2012]{logistic}). More recently, X-ray and optical studies
of some of the brightest of these have identified some strong
``black widow'' candidates: eclipsing binary MSPs with very compact, almost
circular orbits, where the pulsar is destroying its low-mass companion (\cite[Romani \& Shaw 2011, Romani
2012]{Romani11,Romani12}). These studies derive extremely
precise positions (to $\sim0.1^{\prime\prime}$) and stringent constraints on two
out of the three requisite orbital parameters (leaving only the projected semi-axis relatively
unconstrainted). Thus, for the first time, searches for such
binary MSPs using $\gamma$-ray data are possible. A deep LAT search
of 2FGL J2339.6+0532 (the most promising of these ``black
widow'' candidates) unfortunately produced no significant candidate,
but further efforts are in progress (see Belfiore (2012) for details).

\section{Searches for pulsars around the Galactic Center}

Understanding the $\gamma$-ray emission from the Galactic Center (GC)
region is both challenging and controversial. Claims that $\gamma$-ray emission from the GC may be related to dark
matter (e.g. \cite[Hooper \& Goodenough 2011, Weniger 2012]{HG11,Weniger12}) should be measured up against more conventional explanations, such as a possible origin from $\gamma$-ray pulsars.
As a site of massive star formation, it likely contains thousands of
pulsars, but their radio detection is hampered by the large
amount of interstellar scattering. Nevertheless, some pulsars have been discovered fairly close to
the GC, and predictions for the number of radio pulsars that could be
associated with the GC range in the thousands (\cite[Deneva et al. 2009]{Deneva09}). It
is possible that some of the undiscovered pulsars in the
region are {\it radio quiet}, like the majority of young pulsars found in
LAT blind searches (e.g. PSR J1732-3131, see Figure 3).

Blind searches for $\gamma$-ray pulsars around the GC are affected
by low fluxes (due to the large distance), and high levels of diffuse
emission. Indeed, most LAT pulsars are likely
nearby ($\sim$few kpc). The LAT has, however, detected pulsations
from MSP J1823-3021A, in the globular cluster NGC 6624, at 8.4 kpc
(roughly the distance to the GC). It is possible to estimate how far blind searches might
be sensitive out to. From scaling arguments, considering that
the $\gamma$-ray flux of the Crab is several thousand times brighter than the faintest
$\gamma$-ray pulsar discovered in a blind search, we
conclude that it is possible to discover a {\it Crab-like}
$\gamma$-ray pulsar in a LAT blind search out to at least $\sim$15 kpc.

Perhaps the biggest complication in searching for faint pulsars around
the GC is the fact that such a pulsar may be young and {\it noisy} (or
``glitchy''). As described above, the loss of coherence of the signal limits the amount of data
that can be effectively searched (see Figure 2). Thus, regardless of
the number of years of LAT data available, it may only be possible to
search several months at a time. Thus, the sensitivity of searches in this region may only improve significantly
with improvements in reconstruction (e.g. Pass 8, see Baldini (2011),
{\it Fermi} Symposium), or a change in the observing mode. For the first
four years of its mission, {\it Fermi} has operated mostly in survey
mode. This has many advantages for most of the scientific goals of the mission, including pulsar searches (and
timing). Modifying the observing profile, however, could enhance the
sensitivity to detection of faint signals (both from pulsars and/or dark
matter) from the GC. Figure 3 (right) shows the relative gain in
exposure (compared to survey mode) with a modified mode in which the
LAT points to a location slightly offset from the GC. Whenever the GC is
occulted, the LAT could go back into survey mode, thus maintaining
some level of exposure over the entire sky.

\begin{figure}[htbp]
\includegraphics[height=.19\textheight]{GC1.eps}
\includegraphics[height=.20\textheight]{gainmap_75.eps}
\caption{{\bf Left} LAT view of the GC (Porter, 2011 {\it Fermi} Symposium). The bright source at the GC (2FGL J1745.6-2858) is curved and
  non-variable, like most LAT pulsars. {\bf Right}
  Relative exposure (compared to survey mode) with a LAT viewing
  strategy favoring the GC. The gain in exposure is a factor of
  $\sim$3, while the exposure for other parts of the sky decreases to
  roughly one third, in the worst case. }
\label{fig6}
\end{figure}

\section{Conclusions and Prospects}
The NASA Senior Review recently recommended that {\it Fermi} operations
continue through 2016. In principle, {\it Fermi} could continue to
operate well beyond this date.\footnote{Unlike its predecessor,
  EGRET, the LAT has no consumables that limit its lifetime.} As the
LAT accumulates more data, it will detect $\gamma$-ray
pulsations from ever fainter {\it radio-loud} pulsars. As for blind
searches of $\gamma$-ray data, they too will continue to produce new
discoveries, although the increasing computational demands will
require additional resources (e.g. Einstein@Home\footnote{\tt{http://einstein.phys.uwm.edu}}) and
efficient computational techniques to exploit them.
Finally, multi-wavelength observations (radio, X-rays, and optical),
will continue to play a crucial role in LAT pulsar studies, especially
in facilitating the search for more exotic pulsar systems, such as
``black widow'' systems, young binary pulsars, and pulsars around the
GC. In this last case, a modified observing profile
enhancing the exposure to this region (by a factor of $\sim$3) could
play a crucial role in future discoveries.

\section*{Acknowledgements}
The $Fermi$ LAT Collaboration acknowledges support from a number of
agencies and institutes for both development and the operation of the
LAT as well as scientific data analysis. These include NASA and DOE in
the United States, CEA/Irfu and IN2P3/CNRS in France, ASI and INFN in
Italy, MEXT, KEK, and JAXA in Japan, and the K.~A.~Wallenberg
Foundation, the Swedish Research Council and the National Space Board
in Sweden. Additional support from INAF in Italy and CNES in France
for science analysis during the operations phase is also gratefully
acknowledged. I acknowledge the support of the American Astronomical
Society and the National Science Foundation in the form of an
International Travel Grant, which enabled me to attend this conference.

\def \apjl {ApJL}
\def \apj {ApJ}
\def \apjs {ApJS}
\def \jcap {JCAP}