Fukui matrix
The Fukui function is defined as \[\begin{equation} \require{physics} f( \vb{r} ) % = \qty( \pdv{ \rho( \vb{r} ) }{N} )_{ v( \vb{r} ) } \end{equation}\] where \(\rho( \vb{r} )\) is the electron density which depends on the coordinates but implicitly on the number of electrons \(N\) too. Reminiscent of the relation between the \(N\)-electron probability density and density matrix, one might generalize the Fukui function by introducing a Fukui matrix (Bultinck et al. 2011): \[\begin{equation} \require{physics} f( \vb{r}, \vb{r'} ) = \qty( \pdv{ \rho( \vb{r}, \vb{r'} ) }{N} )_{v(\vb{r})=v(\vb{r'})} \end{equation}\] where the requirement \(v(\vb{r})=v(\vb{r'})\) rises from the limitation of constant external potential. In basis-set representation, the density matrix can be rewritten as \[\begin{equation} P_{pq} = \matrixel{\phi_p(\vb{r})}{\rho(\vb{r}, \vb{r'})}{\phi_q(\vb{r'})} \end{equation}\] such that the Fukui matrix projected on a given basis set becomes \[\begin{equation} \vb{f} = \qty{ f_{pq} = \qty( \pdv{P_{pq}}{N} )_{v(\vb{r})=v(\vb{r'})} } \thinspace . \end{equation}\] Since the left and right limit of the Fukui functions are defined as differences of two density matrices with different number of electrons \[\begin{align} \vb{f}^+ &= \vb{P}^{N+1} - \vb{P}^N \\ \vb{f}^- &= \vb{P}^N - \vb{P}^{N-1} \thinspace , \end{align}\] we need to express both density matrices in the same orbital basis. At the Kohn-Sham level, all matrices can be expressed in terms of the real MOs of the neutral molecule. The transformation of the ionic species requires a straightforward unitary transformation \[\begin{equation} \phi^{\text{ion}}_p = \sum_q^K \phi^{\text{neutral}}_q T_{qp} = \sum_q^K \phi^{\text{neutral}}_q (C^{\text{neutral}}_{qp})^{-1} C^{\text{ion}}_{qp} \end{equation}\] and thus the transformation matrix is given by \((\vb{D}^{-1} \vb{C})\) \[\begin{align} \vb{f}^+ %&= \vb{C}^\top \vb{D} \vb{P}^{N+1} \vb{D}^\top \vb{C} %- \vb{P}^N %\\ &= (\vb{C}^{-1} \vb{D}) \vb{P}^{N+1} (\vb{C}^{-1} \vb{D})^\top - \vb{P}^N \\ \vb{f}^- %&= \vb{P}^N %- \vb{C}^\top \vb{D} \vb{P}^{N-1} \vb{D}^\top \vb{C} %\\ &= \vb{P}^N - (\vb{C}^{-1} \vb{D}) \vb{P}^{N-1} (\vb{C}^{-1} \vb{D})^\top \end{align}\] where we have opted for a simpler notation and renamed \(C^{\text{ion}}\) and \(C^{\text{neutral}}\) to \(\vb{D}\) for the coefficient matrix of the ionic species and \(\vb{C}\) for that of the neutral species respectively. The density matrix \(\vb{P}\) in the equations above is always expressed in terms of its own molecular orbitals. Owing to the normalization of the Fukui function, the trace must be equal to 1. Diagonalization of the Fukui matrix results in the Fukui eigenvalues and its associated Fukui orbitals or Fukui naturals.
Physical interpretation
At the Kohn-Sham level, there is always one dominant eigenvector with an eigenvalue equal to 1. The remaining eigenvalues are either zero or come in pairs with opposite sign. Analysis of the Frontier Molecular Orbital (FMO) coefficient gives information of the quality of the FMO approximation. If the coefficient is not close to 1, the Fukui orbital has more contributions than only the FMO approximation and further analysis based on the FMO approximation is ill-advised. The Fukui matrix also tells us whether a negative Fukui function will appear in a certain region, based on negative eigenvalues.