I think this is the same answer as yasmar's, but from a more computational perspective.
It makes sense that the pair of vectors in the two subspaces forming the least angle with each other is distinguished and these vectors should be rotated into each other if the interpolation is to follow a geodesic. To find these vectors, you can find the stationary values of the dot product between arbitrary vectors in the subspaces,
$$(\cos\phi\vec{a}_1+\sin\phi\vec{a}_2)\cdot(\cos\theta\vec{b}_1+\sin\theta\vec{b}_2)\;,$$
where $\vec{a}_1$ and $\vec{a}_2$ form an orthonormal basis for the first subspace and likewise $\vec{b}_1$ and $\vec{b}_2$ for the second one. (You can easily calculate these by orthonormalization if you don't already have them.) With $g_{ij}:=\vec{a}_i\cdot\vec{b}_j$, this becomes
$$
g_{11}\cos\phi\cos\theta+
g_{12}\cos\phi\sin\theta+
g_{21}\sin\phi\cos\theta+
g_{22}\sin\phi\sin\theta\;.
$$
Differentiating with respect to $\theta$ and $\phi$ yields
$$
\begin{eqnarray}
-g_{11}\cos\phi\sin\theta
+g_{12}\cos\phi\cos\theta
-g_{21}\sin\phi\sin\theta
+g_{22}\sin\phi\cos\theta
&=&0\;,\\
-g_{11}\sin\phi\cos\theta
-g_{12}\sin\phi\sin\theta
+g_{21}\cos\phi\cos\theta
+g_{22}\cos\phi\sin\theta
&=&0\;,
\end{eqnarray}
$$
and adding and subtracting these equations leads to
$$
\begin{eqnarray}
(g_{22}-g_{11})\sin(\phi+\theta)+(g_{12}+g_{21})\cos(\phi+\theta)&=&0\;,\\
(g_{22}+g_{11})\sin(\phi-\theta)+(g_{12}-g_{21})\cos(\phi-\theta)&=&0\;,\\
\end{eqnarray}
$$
from which $\phi+\theta$ and $\phi-\theta$, and hence $\phi$ and $\theta$, can be determined.
That gives you the vectors $\vec{a}_<:=\cos\phi\vec{a}_1+\sin\phi\vec{a}_2$ and $\vec{b}_<:=\cos\theta\vec{b}_1+\sin\theta\vec{b}_2$ that should be rotated into each other. The vectors $\vec{a}_>$ and $\vec{b}_>$ with $\vec{a}_> \perp \vec{a}_<$ and $\vec{b}_> \perp \vec{b}_<$ also fulfill $\vec{a}_> \perp \vec{b}_<$ and $\vec{b}_> \perp \vec{a}_<$ (else the stationary point wouldn't be stationary). You can choose their signs so as to form the smaller of the two possible angles between them.
Now these pairs of vectors form two planes that are orthogonal to each other and within which the interpolation rotations should take place, with $\vec{a}_<$ being rotated into $\vec{b}_<$ and $\vec{a}_>$ being rotated into $\vec{b}_>$, both with (different) linear rates of change of the respective rotation angles.
Best Answer
Suggestion: Move your origin to $p$, orthonormalize $m,n$ to get $e_1$,$e_2$ and supplement with two orthonormal vectors $e_3,e_4$. Let $E=[e_1\ e_2\ e_3 \ e_4]$ denote the (4 by 4) matrix of these column vectors.
A vector in these coordinates is described by $v=Ex$ with $x\in {\Bbb R}^4$.
A rotation in the plane spanned by $e_1,e_2$ is given by $f_a(Ex)=E M_a x=E M_a E^{-1} v= E M_a E^T v$ with $$ M_a=\left( \begin{matrix} \cos a & -\sin a & 0 & 0\\ \sin a & \cos a & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{matrix} \right)$$