First, we rewrite the sum as a sum over the full period
$$
S(a,N)=\frac{2}{N+1}\sum_{j=1}^{N+1} \sin^2\left( \frac{2\pi j}{N+1} \right)\sin \left( 2 a \cos \left( \frac{2\pi j}{N+1} \right) \right).
$$
Denote
$$
g(k,m)=\sum _{j=1}^m \cos ^{k}\left(\frac{2 \pi j}{m}\right).
$$
Then the coefficients $a^{2k+1}$ of the Taylor series expansion of $S(a,N)$ is
$$
\frac{2(-1)^k}{N+1}\frac{2^{2k+1}}{(2k+1)!}\left(g(2k+1,N+1)-g(2k+3,N+1)\right).
$$
It is known that (see e.g. this article )
$$
g(k,m)=\frac{m}{2^k}\sum_{r=-\lfloor k/m\rfloor,\,rm+k\,\mathrm{even}}^{\lfloor k/m\rfloor}\binom{k}{\tfrac{rm+k}{2}}.
$$
In particular
$$
g(2k+1,2n+1)=0,\qquad k=0,1,2,\ldots,n-1.\quad (n=N/2)
$$
Observe that this is also a consequence of Chebyshev-Gauss quadrature.
Now, summation with the above formula shows that coefficient of $a^{2k+1}$ in the Taylor series expansion of $S(a,N)$ with $k=0,1,2,\ldots,n-2$ vanish. The first non-vanishing term corresponds to $k=n-1$:
\begin{align}
S(a,N)&\approx (-1)^{N/2}\frac{2(2 a)^{N-1}}{N!}g(N+1,N+1)\\
&=(-1)^{N/2}\frac{2(2 a)^{N-1}}{N!}\frac{2(N+1)}{2^{N+1}}\\
&\approx \frac{(-1)^{N/2} a^{N-1}}{(N-1)!}.
\end{align}
This approximate formula has been confirmed numerically: With increasing $N$, the ratio of the sum and the approximate expression seems to approach $1$.
The higher order coefficients probably can be estimated too with the formulas above, and they will be much smaller than the leading coefficient. This is left as an exercise.
