The clearest explanation of the Clausius-Mossotti (CM) relation I have ever come across is this paper, i.e. Local-field effects and effective-medium theory: a microscopic perspective, D. E. Aspnes, Am. J. Phys. 50, 704 (1982). (I apologise that I can only find a version that is behind a paywall.) The correct definition of the dipole moment must always relate to the microscopic field acting on the individual lattice sites. It is this microscopic field which induces the dipole moments. The microscopic field is different from the apparent macroscopic externally applied field. The latter is the sum of the microscopic applied field and the volume averaged dipole field (related to the macroscopic polarisation field). This is exactly what your second source is saying with $$\mathbf{E}_{eff} = \mathbf{E} + \frac{\mathbf{P}}{3\epsilon_0}.$$ $\mathbf{E}_{eff}$ is the microscopic field acting on each dipole, written in terms of the macroscopically averaged electric field $\mathbf{E}$ and polarisation $\mathbf{P}$. The factor $\frac{1}{3}$ accompanying $\mathbf{P}$ arises due to the volume averaging. I don't know the details of Griffith's derivation, but his symbol $\mathbf{E}$ must denote this microscopic field also, or he has done something dodgy.
The rest of your confusion appears to stem from definitions and units. You are free to define the polarisability $\alpha$ so that $\mathbf{p} = \alpha \epsilon_0 \mathbf{E}$ or so that $\mathbf{p} = \alpha^{\prime} \mathbf{E}$. You convert from one to the other by $\alpha^{\prime} = \alpha \epsilon_0$, exactly as you convert between your corresponding CM expressions. The appearance of a $\frac{1}{4\pi}$ in place of $\epsilon_0$ is common when converting from SI units to other unit systems common in electromagnetism, where often $\epsilon_0 = 1$ by definition.
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
I think I can clear up most of this. Maybe someone whose qm chops are better than mine could help with the parts I'm fuzzy on.
Suppose an even-even nucleus has a prolate deformation (like an American football). This is very common, and basically occurs for any nucleus whose N and Z are both far from any magic numbers. What we really mean when we say that it's deformed is that there are correlations between the different nucleons (neutrons and protons), and the correlations have a certain spatial pattern or organization.
But in its ground state, this nucleus has spin 0. A spin $I=0$ has only a single $I_z$ state, which is $I_z=0$. That means that a zero spin has no orientation degree of freedom. Now how can this be when the thing is supposed to be shaped like a football? Obviously it's possible to give a football different orientations, and those orientations are distinguishable. Well, the general idea is to think of the spin-0 ground state as a superposition of every possible orientation. (I don't know if this description is really rigorous, but I think it's good enough for the present purpose. We have issues like these, and also states with similar orientations can have nonvanishing inner products with each other.)
So in the ground state, the neutrons and protons have this quadrupole-type pattern of correlations with each other, but they have no such correlation with anything external.
I think the way this shows up in the formula you posted is that for $I=0$, it misbehaves. (The numerator and denominator are both zero.)
The formula also misbehaves if you plug in, e.g., $I=1/2$, $I_z=1/2$, $\eta=0$, and $I^2\rightarrow I(I+1)=3/4 $. Here I don't know if there's any geometrical explanation as simple as the one I gave above. This may be where you can't avoid the Wigner-Eckart theorem.
Let's come back to the deformed even-even nucleus with a spin 0 ground state. Although you can't orient the ground state, you can take this nuclear shape and excite it into a state of end-over-end rotation. If you do this, you get a rotational band with spins 0, 2, 4, ... and energies that go approximately like $E\propto I(I+1)$, which is basically the classical $I^2$ result for a rotor, with a quantum correction term added on. The existence of a set of states with this pattern of spins and energies is one of the classic ways of verifying that the nucleus really is deformed. (A spherical nucleus can't rotate collectively.) But in NMR or NQR you would never see the excited states.
In such a rotational band, we also observe anomalously fast electromagnetic transitions such as $4\rightarrow2$ and $2\rightarrow0$. These transitions are fast because they arise from the collective rotation of the whole nucleus, which makes it radiate coherently like an antenna. The speed of these transitions can be described using a transition quadrupole moment, which is different from, but related to, the static quadrupole moment of the nuclear shape. In the excited states, the nucleons have quadrupole correlations not just with each other, but also with the outside world. We can see this because the gamma radiation emitted in the transitions is externally detectable and also has a radiation pattern that is asymmetric as measured in the lab.
When nuclear physicists say that the spin 0 ground state of a rotational band "has" a certain quadrupole moment, what we really mean is that we make a somewhat unrealistic model in which we break rotational symmetry (and therefore slightly violate angular momentum conservation) by modeling the nucleus as if it had a fixed orientation in the lab frame. In this model, the nucleus has a quadrupole moment.