I am interested in decomposing an arbitrary rotation in 3D space into the product of two rotations which are constrained to have their axes in the same plane (for instance x-y plane).
Statement of the conjecture
More precisely stated: Let $R$ be an arbitrary 3D rotation. We want to prove that $R_1$ and $R_2$ exist such that $R=R_1R_2$ with $R_1$ and $R_2$ being rotations with their axes in the, say, $x-y$ plane. The rotation angle of $R_1$ and $R_2$ can be arbitrarily chosen.
Empirical evidence
I attempted to verify the problem numerically. I have drawn random numbers which specify the rotation matrices $R_1$ and $R_2$. I then represented on a 3D sphere the rotation obtained by the product $R_1R_2$. The distance of the points represents the rotation angle. The direction identified by the line connecting the point and the origin denotes the orientation of the rotation axis. One sees that the volume is entirely filled (note that only the left half volume is displayed for the sake of clarity, and that the colors encode the distance from the origin).
This supports the conjecture that any arbitrary 3D rotation can be reached with the product of $R_1 R_2$.
Additional remarks
This problem has a direct connection to a physical system. A positive answer to this question would mean that combining together two birefringent plates allows one to realize any arbitrary transformation of the light polarization by controlling the orientation of the plates and the amount of birefringence for each of them.
Note also that it is well clear to me that any 3D rotation can be described with three parameters, for instance, the Euler angles. Here we have $4$ parameters in total which specify the product $R_1 R_2$. Therefore at least one parameter cannot be free.
It is clear to me as well that an arbitrary rotation $R$ can be decomposed as the product $R_x(\phi_1)R_z(\phi_2)R_x(\phi_3)$ where $R_i(\phi)$ is a rotation with respect to $i$-axis by angle $\phi$. Instead, the question posted here is about to decompose a $R$ with two rotations only.
Best Answer
Your given rotation $R(l,\theta)$ around line $l$ by angle $\theta$ is a composition of two symmetries wrt planes $p_1$ and $p_2$ which intersect along $l=p_1\cap p_2$ at angle $\theta/2$:
$$ R(l,\theta) = S(p_1)\circ S(p_2) = S(p_1)\circ S(p)\circ S(p)\circ S(p_2) $$
where $p$ is the plane you are interested in (e.g., $xy$) because $S\circ S=\text{Id}$.
Now, let $l_i=p_i\cap p$ and $\theta_i=2\times\angle(p_i,p)$ be double the angle between $p_i$ and $p$ ($i=1,2$). (You can always select $p_i$ so that neither is parallel to $p$).
Then $$ \begin{align} S(p_1)\circ S(p) &= R(l_1,\theta_1) \\ S(p)\circ S(p_2) &= R(l_2,-\theta_2) \end{align} $$ and
$$R(l,\theta)=R(l_1,\theta_1)\circ R(l_2,-\theta_2)$$