I am reading the Doctorate thesis 'Zero-Field Anisotropic Spin Hamiltonians in First-Row Transition Metal Complexes: Theory, Models and Applications' (link: https://tel.archives-ouvertes.fr/tel-00608878/document).
On page 34 he writes the zero-field splitting Hamiltonian: $$\hat{H}=\hat{\mathbf{S}} \cdot \bar{D} \cdot \hat{\mathbf{S}}=\left(\begin{array}{ccc} \hat{S}_{x} & \hat{S}_{y} & \hat{S}_{z} \end{array}\right) \cdot\left(\begin{array}{ccc} D_{x x} & D_{x y} & D_{x z} \\ D_{x y} & D_{y y} & D_{y z} \\ D_{x z} & D_{y z} & D_{z z} \end{array}\right) \cdot\left(\begin{array}{c} \hat{S}_{x} \\ \hat{S}_{y} \\ \hat{S}_{z} \end{array}\right)$$ He then states: "and the $\hat{S}_x$ and $\hat{S}_y$ operators are replaced by the adequate linear combinations of the $\hat{S}_+$ and $\hat{S}_-$ operators. By applying this Hamiltonian on the basis of the model space (in this case all the $|S,M_S\rangle$ components of the ground state) the interaction matrix is constructed."
On page 49 he gives the following example for a spin-1 state, i.e. for $|S,M_S\rangle=|1,M_S\rangle$: $$\begin{array}{r|ccc} \hat{H}_{\text {mod}} & |1,-1\rangle & |1,0\rangle & |1,1\rangle \\ \hline\langle 1,-1| & \frac{1}{2}\left(D_{x x}+D_{y y}\right)+D_{z z} & -\frac{\sqrt{2}}{2}\left(D_{x z}+i D_{y z}\right) & \frac{1}{2}\left(D_{x x}-D_{y y}+2 i D_{x y}\right) \\ \langle 1,0| & -\frac{\sqrt{2}}{2}\left(D_{x z}-i D_{y z}\right) & D_{x x}+D_{y y} & \frac{\sqrt{2}}{2}\left(D_{x z}+i D_{y z}\right) \\ \langle 1,1| & \frac{1}{2}\left(D_{x x}-D_{y y}-2 i D_{x y}\right) & \frac{\sqrt{2}}{2}\left(D_{x z}-i D_{y z}\right) & \frac{1}{2}\left(D_{x x}+D_{y y}\right)+D_{z z} \end{array}$$
Now my question is, how do the matrix elements become a linear combination of the parameters $D_{ij}$? If someone understands it, could you calculate a single matrix element to show how it works?