Suppose we have $\mathbf{AB} = \mathbf{BA}$, where $\mathbf{A},\mathbf{B} \in SO(3)$.
What facts does this imply about $\mathbf{A}$ and $\mathbf{B}$? Clearly $\mathbf{A} = \mathbf{I}_3$ and $\mathbf{B} = \mathbf{I}_3$ are always solutions.
Given $\mathbf{A} \neq \mathbf{I}_3$, what are the resulting constraints on $\mathbf{B}$? There are certain cases such as where both have rotation angle $\pi$ and perpendicular, coincident, or negative axes of rotation (Orthogonal axes of rotation of angle $\pi$ implies AB=BA), but what are the general constraints?
We may show for example that for any $m \times m \ \textbf{square}$ matrices $\mathbf{X,Y}$, the equation $\mathbf{XY} = \mathbf{YX}$ implies;
$$ (\mathbf{Y}^{\top}\otimes\mathbf{I}_m - \mathbf{I}_m\otimes\mathbf{Y})\mathbf{X}^s = \mathbf{0}$$
where $\otimes$ is the Kronecker product and $\{\cdot\}^s$ is the column stacking operator. This gives us a statement about $\mathbf{X}$ in terms of the nullspace of an $m^2 \times m^2$ matrix with coefficients dependent on $\mathbf{Y}$, but I am so far unable to constrain it further to $SO(3)$.
Let $\mathbf{a} \in \mathbb{R}^3$ be a vector colinear with the axis of rotation of $\mathbf{A}$. Then, $\mathbf{ABa} = \mathbf{Ba}$.
Then $\mathbf{Ba}$ must also be co-linear with the rotation axis of $\mathbf{A}$, $\mathbf{Ba} = k\mathbf{a}, k \in \mathbb{R}$.
Since $\mathbf{B}$ is a rotation matrix and hence scale preserving, $|k| = 1$. This has only two solutions, $k = \pm 1$. Hence $\mathbf{B}$ can correspond to identity, or rotation by $\pi$ about any axis orthogonal to $\mathbf{a}$. By symmetry, $\mathbf{A}$ must also rotate by $\pi$.