The theory of the spin–lattice relaxation of dilute concentrations of paramagnetic ions (having non-degenerate energy levels) in both harmonic and anharmonic insulating crystals is formulated in terms of Green functions, and the Abrikosov diagram technique for spins is used to analyze these Green functions. A set of rate equations for the populations of the spin energy levels of a spin subjected to a low-frequency longitudinal, external magnetic field is derived using intuitive arguments, and the result of the Green function analysis is shown to be in agreement with these equations. Expressions for the transition probabilities are derived which are correct to all orders in the spin–phonon coupling, and to second order in the three-phonon anharmonic coupling. The transition probabilities for the direct process, the Orbach process, the first- and second-order Raman processes, and the anharmonic Raman process are among those obtained. The perturbation series for the transition probability is shown to be an expansion in powers of the parameters [Formula: see text] and [Formula: see text]where f′ and f″ are spin–phonon coupling constants (having dimensions of energy), εD is the Debye energy, [Formula: see text] is the mean square atomic displacement, and a0 is the lattice constant. The diagrammatic method used here has three major advantages over the decoupling procedures previously used to analyze spin–phonon Green functions. (1) Approximations are made by expanding in powers of a small parameter, (2) contact is made with an intuitive formulation of the rate equations, and (3) formulas valid to all orders in the spin–phonon interaction are obtained. An appendix is added which clarifies some existing ideas concerning the Orbach process.