Tight binding summary
When computing the electronic bands of a crystal in the tight-Binding approximation, the standard way to do it is to construct Bloch solution as
$$ \Psi_n(\textbf{r}, \textbf{k}) = \sum_{\textbf{R}}e^{i\textbf{k}\cdot \textbf{R}}\psi_n(\textbf{r}-\textbf{R}) $$
where in general $\psi_n(\textbf{r}-\textbf{R})$ is a LCAO, namely
$$ \psi_n(\textbf{r}-\textbf{R}) = \sum_{n'= 1}^N C_{nn'} \phi_{n'}(\textbf{r}-\textbf{R})\;\;\; n = 1, ..., N $$
and $\phi_n(\vec{r})$ is the $n^{th}$ atomic orbital. To compute the bands, we use
$$ E_i(\textbf{k}) = \frac{\langle\Psi_i(\textbf{r}, \textbf{k})|\hat{H}|\Psi_i(\textbf{r}, \textbf{k})\rangle}{\langle\Psi_i(\textbf{r}, \textbf{k})|\Psi_i(\textbf{r}, \textbf{k})\rangle} $$
if we now introduce the functions
$$ \Phi_n(\textbf{r}, \textbf{k}) = \frac{1}{\sqrt{N}}\sum_{\textbf{R}}e^{i\textbf{k}\textbf{R}}\phi_{n}(\textbf{r}-\textbf{R}) $$
we can define the matrices $H_{ij} = \langle\Phi_i|\hat{H}|\Phi_j\rangle$ and $S_{ij} = \langle\Phi_i|\Phi_j\rangle$ and the eigenvalue equation gives the secular system
$$ E_i = \frac{\sum_{n, n'}H_{nn'}C^*_{in}C_{in'}}{\sum_{n, n'}S_{nn'}C^*_{in}C_{in'}} $$
by minimizing the $E_i$ we an find the eigenvalues, namely the electronic bands.
My questions
Why are we taking the number of LCAO combinations the same as the number of atomic orbitals we are taking into account? Is it because we are looking for $N$ linearly independent functions, therefore we can't make more than that, having $N$ linearly independent atomic orbitals?
If $\phi_n(\textbf{r})$ are the atomic orbitals, namely eigenvectors of the atomic Hamiltonian, they are an orthogonal basis. Then I expect the functions $\Phi_n(\textbf{r}, \textbf{k})$ to be orthogonal as well. If this is true, why do we bother computing $S_{ij}$ when it is just $\delta_{ij}$?