Choosing the rotation angle-rotation axis parameters instead of the conventional Euler angle parameters {α, β, γ} to parametrise rotations, and a spherical tensor operator basis, Happer was the first to propose a novel multipole operator expansion of the rotation operator [W. Happer, Ann. Phys. (N.Y.) 48, 579 (1968)]. As Happer pointed out, such an expansion in irreducible spherical tensor operators readily lends itself to the methods of Racah algebra, a feature we exploit to simplify the calculus of rotation matrix products. Working with a Clebsch–Gordan coefficient expansion [W. Happer, ibid. cit.; M.S. Marinov, Yad. Fiz. 5, 943 (1967)] of the corresponding rotation matrices , and closure relations for these (2J + 1)-dimensional irreducible representations of the rotation group, we state, and prove for the first time, simple recoupling coefficient expansions for rotation matrix products. With an eye towards applications in NMR, where multiple rotations in spin or coordinate space often require the efficient evaluation of such rotation matrix products, we focus on expansions defining the matrices representing the result of two or three consecutive rotations. We show that the coefficients of these expansions are defined by 6 − j symbols and bipolar spherical harmonics (two rotations), or by 9 − j symbols and tripolar spherical harmonics (three rotations). We show that these recoupling coefficient expansions, in conjunction with Racah algebra, are the cornerstone of an unusually simple and direct method of (1) proving the Euler–Rodrigues (quaternion) parameter composition rules for two and three consecutive rotations and (2) constructing rotational invariants expressed in Weyl-invariant form. These methods are distinguished by their independence on the usual algebraic approaches, either those of operator algebra, or ordinary algebra.