In this work we present a coupled cluster based approach to the computation of the spin orbit coupling matrix elements. The working expressions are derived from the quadratic response function with the coupled cluster parametrization, using the auxiliary excitation operator S. Systematic approximations are proposed with the CCSD and CC3 levels of theory. The new method is tested by computing lifetimes of several electronic states of Ca, Sr and Ba atoms, with Gaussian and Slater basis sets. The results are compared with available theoretical and experimental data.