Let me remind you that the Fock space of multiple fermions is defined to be the antisymmetric (fermionic) subspace of the full Fock space
$$
\Gamma_a=\bigoplus_{n=0}^{\infty}H^{\wedge n},
$$
where $\wedge$ stands for the antisymmetric tensor product
$$
v_1\wedge\ldots\wedge v_n=\frac{1}{n!}\sum_{p\in P_n}\sigma_p v_{p(1)}\wedge\ldots\wedge v_{p(n)}.
$$
Here $\sigma_p$ is the sign of the permutation $p$ in the group of permutations $P_n$.
Thus the confusion here comes from the fact that $c_a^{\dagger}c_b^{\dagger}|0\rangle\neq|ab\rangle$ as you seem to state.
Recall that the creation and annihilation operators are defined within the occupation number representation, i.e. $c_a^{\dagger}c_b^{\dagger}|0\rangle=|11\rangle$, where the first number denotes the number of fermions in state $a$ while the second denotes the number of fermions in state $b$. On the other hand, states written in the occupation number representation are defined to be properly antisymmetrized (for fermions) many-body basis states, as forced upon us by the particle indistinguishability. Therefore they are defined within the fermionic Fock space. Any textbook will show so, take a look at the first chapter of Many-Body Quantum Theory in Condensed Matter Physics: An Introduction by Bruus and Flensberg for example. For two fermions described via a single particle basis $\{|a\rangle,|b\rangle\}$ one possible choice is:
$$
|11\rangle=\frac{1}{\sqrt{2}}\left(|ab\rangle-|ba\rangle\right).
$$
Therefore
$$
c_a^{\dagger}c_b^{\dagger}|0\rangle=\frac{1}{\sqrt{2}}\left(c_a^{\dagger}|0\rangle_1\otimes c_b^{\dagger}|0\rangle_2-c_b^{\dagger}|0\rangle_1\otimes c_a^{\dagger}|0\rangle_2\right)
$$
The familiar anticommutativity of these operators is now obvious from this from
$$
c_b^{\dagger}c_a^{\dagger}|0\rangle=\frac{1}{\sqrt{2}}\left(c_b^{\dagger}|0\rangle_1\otimes c_a^{\dagger}|0\rangle_2-c_a^{\dagger}|0\rangle_1\otimes c_b^{\dagger}|0\rangle_2\right)=-c_a^{\dagger}c_b^{\dagger}|0\rangle
$$
In fact one of the great advantages of creation and annihilation operators is that they include the antisymmetry (for fermions) of the wave function implicitly.
Dotting this with $\langle x_1x_2|$ we obtain:
$$
\langle x_1x_2|c_a^{\dagger}c_b^{\dagger}|0\rangle=\frac{1}{\sqrt{2}}\left(\phi_a(x_1)\phi_b(x_2)-\phi_b(x_1)\phi_a(x_2)\right).
$$
Thus there is no inconsistency, both representations show that the particles are entangled.
On the other hand dotting $\langle x_1x_2|$ with $c_a^{\dagger}|0\rangle_1\otimes c_b^{\dagger}|0\rangle_2$ would simply produce
$$
\psi(x_1,x_2)=\phi_a(x_1)\phi_b(x_2)
$$
Therefore there is no inconsistency here also, however, as I said, the important thing to remember is that
$$
c_a^{\dagger}c_b^{\dagger}|0\rangle\neq c_a^{\dagger}|0\rangle_1\otimes c_b^{\dagger}|0\rangle_2
$$
Let me give it a shot:
If I interpret this correctly, $\mathbf{F}$ will be the operator for the full spin of the coupled system, $\mathbf{S}$ will be the operator of the electron spin (usually, one would consider $\mathbf{J}$, the spin containing also spin-orbit coupling, but we are on the S-shell, hence no angular momentum) and $\mathbf{I}$ will be the nuclear spin. Then it should hold that $\mathbf{F}=\mathbf{S}+\mathbf{I}$, right?
First, let's have a look at the hyperfine structure Hamiltonian $\mathbf{H}_{hf}$. By construction of $\mathbf{F}$, the eigenstates of $\mathbf{H}_{hf}$ will be eigenstates of $F^2$ and $F_z$. This is just the same as for angular momentum and electron spin (and we construct $\mathbf{F}$ to have this property - this lets us label the eigenstates by the quantum number corresponding to $\mathbf{F}$). Hence the Hamiltonian must be diagonal in the $|F^2,m_F\rangle$-basis. One can also see that $F^2$ commutes with $I^2$ and $S^2$ (and so does $F_z$), since $\mathbf{F}=\mathbf{I}+\mathbf{S}$.
Now we have a look at $\mathbf{H}_B$, the interaction Hamiltonian with a constant magnetic field. We can see that (up to some prefactor) $\mathbf{H}_B=S_z$. Hence the eigenstates of $\mathbf{H}_B$ must be eigenstates of $S_z$ and thus also of $S^2$ and, since the two operators are independent (they relate to two different types of spins, hence the operators should better commute) also to $I^2$ and $I_z$, if you want.
The crucial problem is that $S_z$ and $F^2$ do not commute. Why?
Well: $\mathbf{F}=\mathbf{I}+\mathbf{S}$ hence $F^2=S^2+I^2+2\mathbf{S}\cdot \mathbf{I}$. Now $S_z$ and $\mathbf{S}$ do not commute, because $S_z$ does not commute with e.g. $S_x$, which is part of $\mathbf{S}$. Since $F^2$ commutes with $\mathbf{H}_{hf}$ and $S_z$ commutes with $\mathbf{H}_B$, but not with $F^2$, we have that $\mathbf{H}_{hf}$ does not commute with $\mathbf{H}_B$. This means that $\mathbf{H}_B$ and $\mathbf{H}_{hf}$ cannot be diagonal in the same basis, hence you need to have off-diagonal elements.
In order to see how the matrix representing $\mathbf{H}_B$ looks like in the $|F^2,m_F\rangle$-basis, you can express the $|m_I,m_S\rangle$-basis (in which $\mathbf{H}_B$ is diagonal) in terms of the other basis. This is exactly what equations (4.21) do. These are obtained by ordinary addition of angular momenta. From there, you can construct the unitary transforming the basis $|m_I,m_S\rangle$ into $|F^2,m_F\rangle$ and $\mathbf{H}_B$ will be the diagonal matrix in the basis $|m_I,m_S\rangle$ conjugated with this unitary.
EDIT:
I'm not quite sure whether I understand correctly what your problem is, but let me elaborate: We want to find the Hamiltonian $\mathbf{H}_B$ in the $|m_Im_S\rangle$ basis. In this basis, it is diagonal, because $\mathbf{H}_B$ is essentially $S_z$ (hence commutes with $S_z$) and it must also commute with $I_z$ since $S_z$ and $I_z$ are independent.
If we order the basis according to $|\frac{1}{2},\frac{1}{2}\rangle,|-\frac{1}{2},-\frac{1}{2}\rangle,|\frac{1}{2},-\frac{1}{2}\rangle,|-\frac{1}{2},\frac{1}{2}\rangle$, then, we can just read off the Hamiltonian: The first and fourth vector are eigenvectors to eigenvalue $\mu B$, the others of $-\mu B$ (by definition of $S_z$, since the second component in $|m_Im_S\rangle=|(SI)I_z,S_z\rangle$ tells us the eigenvalue of $S_z$ that the basis vector corresponds to), i.e.
$$ \mathbf{H}_B=\begin{pmatrix}{} \mu B & 0 & 0 & 0 \\ 0 & -\mu B & 0 & 0 \\ 0 & 0 & -\mu B & 0 \\ 0 & 0 & 0 & \mu B\end{pmatrix}$$
Now, as I said, you just have to change the basis. The matrix transforming the above basis into the new basis is given by eqn. (4.21a-d):
$$U:=\begin{pmatrix}{} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \frac{1}{2} & \frac{1}{2} \\ 0 & 0 & -\frac{1}{2} & \frac{1}{2} \end{pmatrix}$$
where the ordering of the $|Fm_F\rangle$-basis is as for $\mathbf{H}$ in your text.
Now calculate $U\mathbf{H}_B U^{\dagger}$ and that should give you the part of $\mathbf{H}$ coming from $\mathbf{H}_B$ in the $|F,m_F\rangle$-basis and this will be exactly what is written in your book.
EDIT 2: I sort of suspected this, so here is some more linear algebra for the problem. I'll use Dirac notion since I suspect you are more familiar with this:
Now suppose you have given two bases $|e_i\rangle$ and $|f_i\rangle$ and suppose they are orthonormal bases. What we want is a matrix $U$ that transforms one basis into the other (I'll call it $U$, since it'll be a unitary - if the bases are not orthonormal, it'll only be an invertible matrix). So we want a matrix such that
$$ |f_i\rangle:=U|e_i\rangle \qquad \forall i$$
How to construct this matrix? Well, given an equation for $|f_i\rangle$ in terms of the $|e_i\rangle$ will give you the i-th row of the matrix. You can also see the matrix elements in Dirac notation:
$$ \langle e_j|U|e_i\rangle=\langle e_j|f_i\rangle $$
In your case, $|e_i\rangle=|m_Im_S\rangle$ and $|f_i\rangle=|F^2,m_F\rangle$. Hence equation (4.21a) will give you the first row of the matrix (the ordering of the basis vectors $|m_Im_S\rangle$ as I proposed above), (4.21c) the second (notice the basis ordering in the matrix $\mathbf{H}$!) (4.21b) the second and (4.21d) the last row of the matrix. Using the equation for the matrix elements above, you should be able to check that with not too much trouble. You can also easily check that $U$ is indeed a unitary (i.e. $UU^{\dagger}=U^{\dagger}U=\mathbb{1}$.
Then we can calculate the matrix elements:
$$ \langle e_i |\mathbf{H}|e_j\rangle=\langle e_i|U^{\dagger}U\mathbf{H}U^{\dagger}U|e_j\rangle=\langle f_i| U\mathbf{H}U^{\dagger}|f_j\rangle $$, which tells you how the matrix looks like in the other basis.
Best Answer
Intuitively, $\langle i| H^{\prime} | j \rangle$ gives the interaction between states $i$ and $j$. For most physical systems, adding a perturbed potential, is same as introducing this interaction between different energy eigenstates. In other words, you can also say that the diagonal element of full Hamiltonian come only from the original Hamiltonian.