I may be wrong but this is my take on it:

a "free" S=5/2 system of d-electrons means a d^5 configuration: every d-orbital is occupied once, all have their spins parallel, it's a pure spin system, no orbital contribution to the magnetic moment and hence g=g0.

If your system isn't like this, it means that you have some (small but finite) configurational admixture of other spin states (of which S=3/2 will be of lowest energy). This part of the state now does have finite ligand/crystal field splitting and spin orbit coupling (soc), and this combination is what gives rise to a variation in g-factor. In general, the g-factor becomes anisotropic as a result.

I can only guess that one cand find a relationship between ligand/crystal field and what you denote as rhombicity. For d-electron-systems, the ligand field term is usually way strongest, so that spin orbit coupling can be treated as a perturbation. But a perturbation to what? The easiest cases would be those of spherical or at least cubic symmetry. Deviation from those then have also to be treated as perturbation. You have then two perturbation terms of different character, and the way/amount to which soc intervenes in the results depends on the amount of distortion (rhombicity). I have no difficulty in believing that this can only work in the limit of small distortion.

And perturbation theory anyway implies that the departure from the reference (unperturbed) state is small.

Mathematically you may think of this as follows: think of a basis set of states of the mixed configuration, symbolically represented as below

{ ams * |S=5/2, mS=-5/2...5/2> + bms' * |S'=3/2, mS'=-3/2...3/2> }

Zeeman energy, ligand field and spin orbit coupling are the ingredients to the total Hamiltonian the electrons are subject to. These do not commute, in general and therefore have no common set of eigenstates. To obtain the energies of all states (and the g-factor is in the energy differences of the lowest energy states at small fields) you would have to diagonalize the total Hamiltonian (compute eigenvalues and eigenvectors of the Hamiltonian). You therefore have to find the roots of the characteristic polynomial of the matrix representing the Hamiltonian. As soon as there are more than 4 states, there is no analalytical solution for the roots of that polynomial, you have to solve the thing numerically.

In perturbation theory of first order, you just consider the effect of the different terms on the diagonal elements of that matrix (instead of diagonalizing it). "Energies change but the states don't" in this approximation. For this to be correct, the states chosen as the basis in which the Hamiltonian is represented need to be as close as possible to the 'true' ones. You cannot afford large deviations. Getting back to your question: if the reference state for the g-factor formulae is one with zero rhombicity, perturbation theory for the Zeeman energy (-> g-factor) gets out of control if the rhombicity gets too large.

Sorry for the long post. And no, I have no suitable reference at hand to recommend.

Hello, thank you for your answer. I understand the method which you propose. I have done that for high spin Co(II). But in this case I'm only interested in derivation of formulae using a full spin Hamiltonian for S = 5/2. Any high spin system may have a zero-field splitting. In the case of Fe(III) I've read of cases in which the D parameter is ~10 cm^-1, but actually I haven't done a full survey of the values it can have. The ZFS does actually come from spin-orbit plus crystal field (plus possibly other interactions). But once a spin hamiltonian is defined, and if it is enough to describe the magnetic properties (which may not always be the case), one can diagonalize only the spin energy matrices. This is done numerically for Fe(III) (S=5/2), but I wanted to know if there are any approximate analytical formulas, as is the case with Co(II) S = 3/2 (Pilbrow, J. Magn. Reson, 1978, 31, 479-490).

Just with respect to your statement "Any high spin system may have a zero-field splitting":

If the spin corresponds to a half filled shell then the state is spherically symmetric. No single electron crystal field splitting can change this as long as you say "the state is d5" and "d-orbitals are close enough to the truth". There is no orbital moment then, and therefore no spin orbit coupling. (the single electron soc operator has zero trace, and therefore with all orbitals occupied once, there is zero net effect)

And in absence of soc g=g0, as you know.

(Pls. do correct me if I am wrong with this)

I meant "may" have, because matematically the term S*D*S does not vanish for S > 1/2. But of course its value could be zero or very close to zero. As I'm only interested in organometallic/biological compounds, there are always distortions that give a non-negligible ZFS. But for example for some Mn(II) compounds with S = 5/2, the ZFS can be around 500 G (~0.05 cm^-1) or even smaller. But for Fe(III) I think (although I´m not sure) that usually the ZFS is not negligible and even larger or much larger than the Zeeman interaction in normal EPR (X-band). That is the kind of compounds in which I'm currently interested.

Thank you very much!