When $N=1$, a global minimiser can be obtained explicitly. The objective function in this case is equal to $\frac12(\|a\|^2+\|b\|^2)-a^TRQR^Tb$ and the problem boils down to maximising $\langle a,Q'b\rangle$ where $Q'=RQR^T$.
If $a=0$ or $b=0$, every $R$ gives an optimal solution. Suppose $a$ and $b$ are nonzero. Normalise them to unit vectors and let $\phi\in[0,\pi]$ be the angle between them. That is, we assume that $\|a\|=\|b\|=1$ and $\phi=\arccos\langle a,b\rangle\in[0,\pi]$. Since $Q$ is a rotation matrix and $Q'$ is orthogonally similar to $Q$, the two matrices must have the same angles of rotation. Therefore
$$
Q=V\pmatrix{\cos\theta&-\sin\theta&0\\ \sin\theta&\cos\theta&0\\ 0&0&1}V^T
\ \text{ and }\ Q'=U\pmatrix{\cos\theta&-\sin\theta&0\\ \sin\theta&\cos\theta&0\\ 0&0&1}U^T
$$
for some $\theta\in[0,\pi]$ and some rotation matrices $U$ and $V$ whose last columns are the rotation axes of $Q$ and $Q'$ respectively. Once $Q'$ is determined, we have $Q'=RQR^T$ by taking $R=UV^T$.
We now consider three possibilities. First, when $\phi=0$, $a$ is equal to $b$. Therefore $\langle a,Q'b\rangle$ is maximised and its value is $1$ when the rotation axis of $Q'$ is $a$. Hence we can pick a $U\in SO(3)$ whose last column is $a$.
Second, when $\phi\ge\theta>0$, since the great circle distance satisfies the triangle inequality, no rotation of $b$ by the angle $\theta$ about any axis can reduce the angle between $a$ and $b$ to a value smaller than $\phi-\theta$. Hence the maximum possible value of $\langle a,Q'b\rangle$ is $\cos(\phi-\theta)$ and this is attained when both $a$ and $b$ lie on the plane of rotation and $Q'$ rotates $b$ towards $a$. So, we can pick a rotation matrix $U$ whose last column is $\frac{b\times a}{\|b\times a\|}$ if $b\ne-a$, or any unit vector orthogonal to $a$ if $b=-a$.
Finally, when $0<\phi<\theta$, by picking an appropriate axis, one can rotate the vector $b$ by an angle $\theta$ to align with $a$. Therefore the maximum value of $\langle a,RQR^Tb\rangle$ is $1$ and it is attained when $RQR^Tb=a$. For example, let $x=\frac{b-a}{\|b-a\|},\, y=\frac{b+a}{\|b+a\|} $ and $z=x\times y$. Then $x,y,z$ form an orthonormal basis of $\mathbb R^3$,
$$
b=\sin\left(\frac{\phi}{2}\right)x+\cos\left(\frac{\phi}{2}\right)y
$$
and the similar holds for $a$, except that the coefficient of $x$ is negated. Now let $y'$ and $z'$ be the results of rotating $y$ and $z$ about the axis $x$ by an angle $\omega$. That is, $y'=\cos(\omega)y+\sin(\omega)z$
$$
z'=\cos(\omega)z-\sin(\omega)y.\tag{1}
$$
Then
$$
b=\sin\left(\frac{\phi}{2}\right)x+\cos\left(\frac{\phi}{2}\right)\cos(\omega)y'-\cos\left(\frac{\phi}{2}\right)\sin(\omega)z'
$$
and the similar holds for $a$ except that the coefficient of $x$ is negated. It follows that $a$ and $b$ have the same components along $z'$ and the angle $\theta'$ between their orthogonal projections on the $xy'$-plane is determined by the equation $\tan\left(\frac{\theta'}{2}\right)=\frac{\tan(\phi/2)}{\cos(\omega)}$. Hence $\theta'=\theta$ if we take
$$
\omega=\arccos\left(\frac{\tan(\phi/2)}{\tan(\theta/2)}\right).
$$
With this $\omega$, we can pick a rotation matrix $U$ whose last column is the $z'$ in $(1)$.
When $N\ge2$ the problem can be reduced to a spherically constrained quadratic program that can be solved iteratively, but I will not go into details here.
