If we express the density matrix in operator form, it is more easy to use in NMR calculation. Let me show you.

for spin ½ ensemble, the density matrix has 4 matrix elements. and we are always to expand it into Pauli matirx and identity matrix, which also give us 4 coefficients.

$\rho = \begin{pmatrix} \rho_+ & \rho_2 \\ \rho_1 & \rho_- \end{pmatrix}$

$\rho = \frac{1}{2} \left( (\rho_++\rho-)1 + (\rho_+-\rho_-) \sigma_z + (\rho_1+\rho_2) \sigma_x + (\rho_1-\rho_2) \sigma_y \right)$

the Pauli matrix is just the angular momentum operator in some sort of reduced form. Thus, the magnetization or the polarization can be very easy to calculated by notice that:

$\sigma_i \sigma_j = i\epsilon_{ijk} \sigma_k$

and the result is the trace.

in a rotation frame, the density matrix under transform:

$\rho_R = R_z(-\omega_0 t) \rho R_z(\omega_0 t)$

and the Pauli matrix change under rotation is like a vector following right-hand rotation rule:

$\sigma_z \rightarrow \sigma_z$

$\sigma_x \rightarrow \sigma_x cos(\omega_0 t) + \sigma_y sin(\omega_0 t)$

$\sigma_y \rightarrow -\sigma_x sin(\omega_0 t) + \sigma_y cos(\omega_0 t)$

Thus, for a polarized ensemble, the off-diagonal matrix element are zero, the density matrix keeps the same under the changing of frame.

in the NMR calculation, the flipping of polarization is just applying the suitable rotation operator. for example, for a $\pi_x$ – pulse, the density matrix becomes:

$\rho_R \rightarrow R_x(-\pi)\rho_R R_x(\pi)$

For the longitudinal relaxation, which is only affect the z-component.

$\sigma_z \rightarrow (1-2 e^{-t/T_1}) \sigma_z$

for the transverse relaxation of the FID ( free induction decay ), which is due to de-synchronization of individual spin ):

$\sigma_x \rightarrow (\sigma_x cos(\omega_0 t) + \sigma_y sin(\omega_0 t))e^{-t/T_2}$

$\sigma_y \rightarrow (-\sigma_x sin(\omega_0 t) + \sigma_y cos(\omega_0 t))e^{-t/T_2}$

since the density matrix is now written on Pauli Matrix, and the polarization can be easily to see :

$P = \left< J \right> = Tr(( J_x , J_y, J_z) \rho) = (\rho_+ - \rho_- , \rho_1+\rho_2 , \rho_1 - \rho_2 )$

thus, we can write it as:

$\rho = f 1 + g \vec{\sigma} \cdot \vec{n}$

with

$f = \frac{1}{2} ( \rho_+ + \rho_- )$

and g is the magnitude of P/4.

by seeing in this way , we can connect it the the matrix M in spin scattering.