Why is C a matrix and not a vector in Roothaan Equations?

I'm trying to figure out why $C$ is a matrix and not a vector in the Roothaan Equations. Below, I have a derivation which I have taken from [this](https://www.colby.edu/chemistry/PChem/notes/Roothaan.pdf) PDF in which, to me, $C$ appears to be a vector. But when I look at other derivations (for example [this derivation](https://nznano.blogspot.com/2018/03/simple-quantum-chemistry-hartree-fock.html)) $C$ is a matrix, although I can't figure out why $C$ is a matrix, I feel it is poorly explained.


So this is my derivation where $C$ is a vector:


Here is the basis function for 1 electron:
$$
\Psi\_a = \sum\_{i=1}^n c\_{ia}\phi\_i
$$


This is the Fock operator for this electron $a$:


$$
f\_a = h\_a + \sum\_{J=1}^{N/2}(2J\_j(a)-K\_j(a))
$$
Above, I'm using $a$ to represent a particular electron, so $J(a)$ and $K(a)$ are, with a loose notation, Coulomb and exchange operators for electron $a$.


So we now have an equation:
$$
f\_a \Psi\_a = \epsilon\_a \Psi\_a \\
f\_a \sum\_{i=1}^n c\_{ia}\phi\_i = \epsilon\_a \sum\_{i=1}^n c\_{ia}\phi\_i \\
$$


This is $1$ equation for $n$ unknown coefficients, so to get $n$ equations we can left multiply each side by each of the $n$ $\phi\_i$ and integrate over all space with respect to electron $a$'s location. This would give us $n$ equations:


$$
\int \phi\_j f\_a \sum\_{i=1}^n c\_{ia}\phi\_i = \int \phi\_j\epsilon\_a \sum\_{i=1}^n c\_{ia}\phi\_i \\
\sum\_{i=1}^n c\_{ia} \int \phi\_j f\_a \phi\_i = \sum\_{i=1}^n c\_{ia} \int \phi\_j\epsilon\_a \phi\_i \\
\sum\_{i=1}^n c\_{ia} F\_{ji} = \epsilon\_a \sum\_{i=1}^n c\_{ia} S\_{ji} 
$$


Where $F$ and $S$ form two $nxn$ matrices, and the matrix equation becomes $FC=SC \epsilon\_a$


$$
F = \begin{bmatrix} 
 F\_{11} & F\_{12} & \dots \\
 \vdots & \ddots & \\
 F\_{n1} & & F\_{nn} 
 \end{bmatrix}
\qquad
S = \begin{bmatrix} 
 S\_{11} & \dots & S\_{1n}\\
 \vdots & \ddots & \vdots\\
 S\_{n1} & \dots & S\_{nn} 
 \end{bmatrix}
C = \begin{bmatrix} 
 C\_{1a}\\
 \vdots\\
 C\_{na} 
 \end{bmatrix}
$$


In this case, $C$ is a vector, and I believe we are solving for only 1 electron $a$ and it's energy $\epsilon\_a$. Keep in mind I believe this is wrong. Above, $\epsilon\_a$ would be a constant. The left side has dimensions $(n\times n) \times (n \times 1) = (n \times 1)$ and the right side of the equation has dimension $(n\times n) \times (n \times 1) \times (1 \times 1) = (n \times 1)$


So if we wanted to solve for two electrons, $a$ and $b$, I believe we would just solve one set of these matrix equations for each electron. The Fock matrix and coefficient matrix would be functions of the electron: $F(a)C(a) = S(a)C(a)\epsilon\_a$ and $F(b)C(b) = S(b)C(b)\epsilon\_b$


Edit: Shown below, I tried to form one matrix equation for the case of two electrons by concatenating the $C$ vector and Fock matrix but I can't make the dimensions work:


$$
F = 
\left(\begin{array}{@{}c|c@{}}
 \begin{matrix}
 F\_{11a} & \dots & F\_{12a} \\
 F\_{n1a} & \dots & F\_{nna}
 \end{matrix}
 & 0 \\
\hline
 0 &
 \begin{matrix}
 F\_{11b} & \dots & F\_{12b} \\
 F\_{n1b} & \dots & F\_{nnb}
 \end{matrix}
\end{array}\right)
C = \begin{bmatrix} 
 C\_{1a}\\
 \vdots\\
 C\_{na}\\
 C\_{1b}\\
 \vdots\\
 C\_{nb}
 \end{bmatrix}
$$


$$FC=SC$$


So in this scenario $C$ is still a vector but I could not figure out how to work in the energy while keeping the dimensions correct because now the energy is a $2 \times 1$ vector $[\epsilon\_a, \epsilon\_b]$.


Anyway I clearly have a major assumption wrong. Could someone please show me how to properly form these equations with $C$ as a matrix?


**EDIT**


Responding to PAEP's answer, while I do know the correct form, I just don't get why that is the form. Thus the need for a derivation.


For example if I multiply out the correct form which you supplied:


$$
FC = \begin{bmatrix} 
 (F\_{1,\bullet} \times C\_{\bullet, 1}) & (F\_{1,\bullet} \times C\_{\bullet, 2}) & \dots \\
 \vdots & \ddots & \\
 (F\_{n,\bullet} \times C\_{\bullet, n}) & & (F\_{n,\bullet} \times C\_{\bullet, n}) 
 \end{bmatrix}
$$


Where $(F\_{i,\bullet} \times C\_{\bullet, i})$ is the dot product of row $i$ of $F$ times column $i$ of $C$


Similarly,


$$
SC = \begin{bmatrix} 
 (S\_{1,\bullet} \times C\_{\bullet, 1}) & (S\_{1,\bullet} \times C\_{\bullet, 2}) & \dots \\
 \vdots & \ddots & \\
 (S\_{N,\bullet} \times C\_{\bullet, n}) & & (S\_{n,\bullet} \times C\_{\bullet, n}) 
 \end{bmatrix}
$$


then,


$$
SC\epsilon = \begin{bmatrix} 
 (S\_{1,\bullet} \times C\_{\bullet, 1})\epsilon\_{11} & & & & \\
 & (S\_{2,\bullet} \times C\_{\bullet, 2})\epsilon\_{22} & & & \\
 & \ddots & & \\
 & & & (S\_{n,\bullet} \times C\_{\bullet, n})\epsilon\_{nn} \\
\end{bmatrix}
$$


Can you derive for me $F\_{i, j}$ in this case? If you look at my original derivation of $F$ it was built from 1 electron only. But this $F$ looks to be built from all electrons yet it has no additional terms.




@Frank, I would try to give you a more detailed answer. This answer is inspired in Lewars book which I strongly recommend you consult.


It is important that you take into account than when you applied the Hartre-Fock self consistent method you are optimizing a single electron spin-orbital but you are considering the rest of the electrons of the atom or molecule you are studying. And as @Martin mentioned, the Fock operators $\hat{F}$ are pseudooperators in the sense that in order to obtain the optimized spin-orbitals you have to modify the orbitals that appear in the operators as part of the optimization process.


Lets start with the Hartree-Fock equations


$$\widehat {\rm{F}}{\kern 1pt} {\psi \_i} = {\varepsilon \_{\rm{i}}}{\kern 1pt} {\psi \_i} $$


Expressing $\psi \_i$ as a linear combination of a set of orbitals


$${\psi \_i} = \sum\limits\_{\nu = 1}^K {{c\_{\nu i}}{\phi \_\nu }}$$


you would obtain


$$\begin{array}{l}\widehat {\rm{F}}{\kern 1pt} \sum\limits\_{\nu = 1}^K {{c\_{\nu i}}{\phi \_\nu }} = {\varepsilon \_{\rm{i}}}{\kern 1pt} \sum\limits\_{\nu = 1}^K {{c\_{\nu i}}{\phi \_\nu }} \\
\sum\limits\_{\nu = 1}^K {{c\_{\nu i}}} \widehat {\mathop{\rm F}\nolimits} (1){\phi \_\nu } = {\varepsilon \_i}\sum\limits\_{\nu = 1}^K {{c\_{\nu i}}} {\kern 1pt} {\kern 1pt} {\phi \_\nu }\end{array}
$$


Considering that you work with K orbitals you will have


$$\begin{array}{c}
\sum\limits\_{\nu = 1}^K {{c\_{\nu 1}}} \widehat {\mathop{\rm F}\nolimits} (1){\phi \_\nu } = {\varepsilon \_1}\sum\limits\_{\nu = 1}^K {{c\_{\nu 1}}} {\kern 1pt} {\kern 1pt} {\phi \_\nu }\\
\sum\limits\_{\nu = 1}^K {{c\_{\nu 2}}} \widehat {\mathop{\rm F}\nolimits} (1){\phi \_\nu } = {\varepsilon \_2}\sum\limits\_{\nu = 1}^K {{c\_{\nu 2}}} {\kern 1pt} {\kern 1pt} {\phi \_\nu }\\
\vdots \\
\sum\limits\_{\nu = 1}^K {{c\_{\nu K}}} \widehat {\mathop{\rm F}\nolimits} (1){\phi \_\nu } = {\varepsilon \_K}\sum\limits\_{\nu = 1}^K {{c\_{\nu K}}} {\kern 1pt} {\kern 1pt} {\phi \_\nu }
\end{array}
$$


where $K \ge N$, K is the number of orbitals in the base and N is the total number of electrons.


Multiplying both members on the left by $\phi\_1^\*$ and integrating over the whole space,


$$
\begin{array}{c}
\sum\limits\_{\nu = 1}^K {{c\_{\nu 1}}} {F\_{1\nu }} = {\varepsilon \_1}\sum\limits\_{\nu = 1}^K {{c\_{\nu 1}}} {\kern 1pt} {\kern 1pt} {S\_{1\nu }}\\
\sum\limits\_{\nu = 1}^K {{c\_{\nu 2}}} {F\_{1\nu }} = {\varepsilon \_2}\sum\limits\_{\nu = 1}^K {{c\_{\nu 2}}} {\kern 1pt} {\kern 1pt} {\kern 1pt} {S\_{1\nu }}\\ 
\vdots \\
\sum\limits\_{\nu = 1}^K {{c\_{\nu K}}} {F\_{1\nu }} = {\varepsilon \_K}\sum\limits\_{\nu = 1}^K {{c\_{\nu K}}} {\kern 1pt} {S\_{1\nu }}
\end{array}
$$


Repiting this operation for $\phi\_2^\*$ , $\ldots$, $\phi\_K^\*$


$$
\begin{array}{c}
\sum\limits\_{\nu = 1}^K {{c\_{\nu 1}}} {F\_{2\nu }} = {\varepsilon \_1}\sum\limits\_{\nu = 1}^K {{c\_{\nu 1}}} {\kern 1pt} {\kern 1pt} {S\_{2\nu }}\\
\sum\limits\_{\nu = 1}^K {{c\_{\nu 2}}} {F\_{2\nu }} = {\varepsilon \_2}\sum\limits\_{\nu = 1}^K {{c\_{\nu 2}}} {\kern 1pt} {\kern 1pt} {\kern 1pt} {S\_{2\nu }}\\
\vdots \\
\sum\limits\_{\nu = 1}^K {{c\_{\nu K}}} {F\_{2\nu }} = {\varepsilon \_K}\sum\limits\_{\nu = 1}^K {{c\_{\nu K}}} {\kern 1pt} {S\_{2\nu }}\end{array}
$$


$$ \ldots $$


$$
\begin{array}{c}
\sum\limits\_{\nu = 1}^K {{c\_{\nu 1}}} {F\_{K\nu }} = {\varepsilon \_1}\sum\limits\_{\nu = 1}^K {{c\_{\nu 1}}} {\kern 1pt} {\kern 1pt} {S\_{K\nu }}\\
\sum\limits\_{\nu = 1}^K {{c\_{\nu 2}}} {F\_{2\nu }} = {\varepsilon \_2}\sum\limits\_{\nu = 1}^K {{c\_{\nu 2}}} {\kern 1pt} {\kern 1pt} {\kern 1pt} {S\_{K\nu }}\\
 \vdots \\
 \sum\limits\_{\nu = 1}^K {{c\_{\nu K}}} {F\_{K\nu }} = {\varepsilon \_K}\sum\limits\_{\nu = 1}^K {{c\_{\nu K}}} {\kern 1pt} {S\_{K\nu }}
\end{array}
$$


where ${F\_{\mu \nu }} = \left\langle {{\phi \_\mu }(1)\left| {\widehat F(1)} \right|} \right.\left. {{\phi \_\nu }(1)} \right\rangle$ and ${S\_{\mu \nu }} = \left\langle {{\phi \_\mu }(1)} \right.\left| {{\phi \_\nu }(1)} \right\rangle$


This operation yields a set of K x K equations known as the Roothaan-Hall equations.


You can write this equation in compact form


$$
 \sum\limits\_{\nu = 1}^K {{c\_{\nu i}}} {F\_{\mu \nu }} = {\varepsilon \_i}\sum\limits\_{\nu = 1}^K {{c\_{\nu i}}} {\kern 1pt} {S\_{\mu \nu }}\quad \mu = 1,2,3, \ldots ,K
$$


where %$\psi\_i (i = 1, 2, \ldots , K)$, so


$${F\_{11}}{\kern 1pt} {c\_{11}} + {F\_{12}}{c\_{21}} + \cdots + {F\_{1K}}{c\_{K1}} = {\varepsilon \_1}\left( {{S\_{11}}{\kern 1pt} {c\_{11}} + {S\_{12}}{\kern 1pt} {c\_{21}} + \cdots {S\_{1K}}{\kern 1pt} {c\_{K1}}} \right) $$
$$ \ldots $$
$${F\_{K1}}{\kern 1pt} {c\_{1K}} + {F\_{K2}}{c\_{2K}} + \cdots + {F\_{KK}}{c\_{KK}} = {\varepsilon \_K}\left( {{S\_{K1}}{\kern 1pt} {c\_{1K}} + {S\_{K2}}{\kern 1pt} {c\_{2K}} + \cdots {S\_{KK}}{\kern 1pt} {c\_{KK}}} \right) $$


In matrix form you will get


$${\mathbf{FC}} = {\mathbf{SCE}}$$


**Addendum: The Fock operator**


The Fock operator $\widehat {\rm{F}}$ that appears in the Hartree-Fock equations


$$\widehat {\rm{F}}{\kern 1pt} {\psi \_i} = {\varepsilon \_{\rm{i}}}{\kern 1pt} {\psi \_i}$$


is a one-electron operator.


This operator can be written as


$$\begin{array}{c}\widehat F(1) = - \frac{1}{2}\nabla \_1^2 - \sum\limits\_\mu {\frac{{{Z\_\mu }}}{{{r\_{\mu 1}}}}} + \sum\limits\_{j = 1}^N {(2{{\widehat {\mathop{\rm J}\nolimits} }\_j} - {{\widehat K}\_j})} \\ = \widehat h(1) + \sum\limits\_{j = 1}^N {(2{{\widehat J}\_j} - {{\widehat K}\_j})} \end{array}$$


where "(1)" just reminds as that the operator is a single electron operator.


The operator $\widehat {\rm{h}}(1)$


$$\widehat {\rm{h}}(1) = - \frac{1}{2}\nabla \_1^2 - \sum\limits\_\mu {\frac{{{Z\_\mu }}}{{{r\_{\mu 1}}}}} $$


is also called the *core-hamitonian* since it represents the energy of an electron that occupies orbital $\phi \_i$ under the attraction of a nuclear *core* (i.e. ignoring the interaction with other electrons).


${\widehat {\mathop{\rm J}\nolimits} \_j}(1)$ corresponds to the *Coulomb operator* that accounts for the classical Coulomb electron-electron respulsion


$${\widehat {\mathop{\rm J}\nolimits} \_j}(1){\kern 1pt} \,{\phi\_i}(1) = {\phi\_i}(1){\kern 1pt} \int {{{\left| {{\phi \_j}(2)} \right|}^2}\frac{1}{{{r\_{12}}}}d{\tau \_2}}$$


whose matrix element $J\_{ij}$ is the *Coulomb integral*


$${J\_{ij}} = \left\langle {{\phi\_i}(1)\left| {{{\widehat J}\_j}(1)} \right|} \right.\left. {{\phi\_i}(1)} \right\rangle$$


$${J\_{ij}} = \int {\phi\_i^\*(1)\,\phi\_j^\*(2)\frac{1}{{{r\_{12}}}}{\phi \_i}(1)\,{\phi \_j}(2)\,d{\tau \_1}d{\tau \_2}} $$


Finally, the *exchange operator* ${\widehat {\mathop{\rm K}\nolimits}\_j}(1)$ is given by


$${\widehat K\_j}\left( 1 \right){\phi\_i}(1) = {\kern 1pt} \,{\phi \_j}(1)\int {\phi \_j^\*(2)\frac{1}{{{r\_{12}}}}{\phi\_i}(2){\kern 1pt} d{\tau \_2}} $$


and its matrix element $K\_{ij}$ is the *exchange integral* (which has no classical counterpart and it is a result that appears when you use a Slater determinant to write the electronic wavefunction)


$${K\_{ij}} = \left\langle {{\phi\_i}(1)} \right.\left| {{{\widehat K}\_j}(1)} \right|\left. {{\phi\_i}(1)} \right\rangle $$


$${K\_{ij}} = \int {\phi\_i^\*(1)\,\phi \_j^\*(2)\frac{1}{{{r\_{12}}}}{\phi \_j}(1)\,{\phi \_i}(2)\,d{\tau \_1}d{\tau \_2}}$$


 4

I think that there is a mistake in the the Roothan-Hall equations.


The matrix equation should read


$$ \mathbf{FC = SC \epsilon}$$
where


$$
\mathbf{F} = \begin{bmatrix} 
 F\_{11} & F\_{12} & \dots \\
 \vdots & \ddots & \\
 F\_{n1} & & F\_{nn} 
 \end{bmatrix}
\quad
\mathbf{S} = \begin{bmatrix} 
 S\_{11} & \dots & S\_{1n}\\
 \vdots & \ddots & \vdots\\
 S\_{n1} & \dots & S\_{nn} 
 \end{bmatrix}
$$
$$
\mathbf{C} = \begin{bmatrix} 
 c\_{11} & \dots & c\_{1n}\\
 \vdots & \ddots & \vdots\\
 c\_{n1} & \dots & c\_{nn} 
 \end{bmatrix}
\quad
\mathbf{\epsilon} = \begin{bmatrix} 
 \epsilon\_{11} & 0 & \dots & 0\\
 0 & \epsilon\_{22} & \dots & 0\\
 \vdots & \vdots &\ddots & \vdots\\
 0 & 0 & \dots & \epsilon\_{nn}\\
 \end{bmatrix}
$$


Probably you can find a detailed outline of the solution of the Roothan-Hall equations in many Quantum Chemistry/Computational Chemistry books. You could try *section 5.2.3. The Hartree-Fock Equations* of E. G. Lewars, Computational Chemistry 2nd ed., Springer (2011).


 1

I think you should try a divide-and-conquer approach. Lets begin with a top-down view to avoid getting tangled up in details.


The Roothaan-Hall equations can be stated in the following form


$$
\mathbf{ FC = SC E }
$$
where $\bf F$ is the Fock-matrix which we obtain by expanding the Fock Operator in a discrete basis set of functions. $ \bf C$ is the coefficient matrix that contains the eigenvector coefficients with respect to this basis. $\bf S$ is the overlap matrix of the basis set and $\bf E$ a diagonal matrix with the eigenvalues.


First, lets assume that our basis is orthogonal such that $\bf S$ can be dropped since it becomes a unit matrix, as this is irrelevant to your question of the dimensionality of $\bf C$.


The form of the equation\* becomes then that of a regular eigenvalue problem
$$
\bf FC = CE
$$
and the question how the dimensionality of F and C are connected is the same question as if applied to any other eigenvalue equation.


Eigenvalue problems are often given with a single column vector, like this
$$
M\vec v = \lambda\_v \vec v
$$
Now if $M$ is quadratic and diagonalizable, with dimension $N \times N$, then we know that we have $N$ eigenvectors. Lets say $N=2$, then we can write the whole system of equations down in a form like this
$$
M(\vec v, \vec w) = (M\vec v, M\vec w) = (\lambda\_v\vec v,\lambda\_w\vec w,)=(\vec v,\vec w)\left( \begin{matrix}\lambda\_v & 0 \\ 0 & \lambda\_w \end{matrix}\right) =\left( \begin{matrix}v\_1 & w\_1 \\ v\_2 & w\_2 \end{matrix}\right)\left( \begin{matrix}\lambda\_v & 0 \\ 0 & \lambda\_w \end{matrix}\right)
$$


It is clear that we can denote the whole system of equations using a coefficient matrix, instead of using a label for each individual eigenvector.
$$\begin{align}
MC &= C\Lambda \\
(M \left(\begin{matrix}C\_{11} \\ C\_{21} \end{matrix}\right), M \left(\begin{matrix}C\_{12} \\ C\_{22} \end{matrix}\right))&= (\lambda\_1 \left(\begin{matrix}C\_{11} \\ C\_{21} \end{matrix}\right), \ \lambda\_2 \left(\begin{matrix}C\_{12} \\ C\_{22} \end{matrix}\right))\\
&=\left( \begin{matrix}C\_{11} & C\_{12} \\ C\_{21} & C\_{22} \end{matrix}\right)\left( \begin{matrix}\lambda\_{1} & 0 \\ 0 & \lambda\_{2} \end{matrix}\right)={C\Lambda}
\end{align}$$


Note that all matrices involved have the dimension of $N\times N$. The overlap matrix would also have the same dimension. You get $N$ eigenvectors and each eigenvector has $N$ coefficients. "Putting" does eigenvectors side-by-side into a matrix yields the coefficient matrix with dimension $N \times N$. The details how the matrix elements are obtained are irrelevant to any of these points.


\* It only looks like a linear eigenvalue equation, but the problem is nonlinear since the matrix elements of the Fock operator are functions of the coefficient matrix elements, which is why we have to solve the problem iteratively until self-consistency is reached.


 1