This is only a very partial answer. Suppose that $a<-n$ and $p$ is a non-integer in the interval $(-1,n)$. Then $2j-n-a>0$ for all $j=0,\dots,n$. So, by the mean-value theorem applied (say, repeatedly) to the $n$-fold symmetric difference in the expression of $I_{n,p,a}$ in $(**)$ in the question statement, one has $$I_{n,p,a}=i^p\, \Gamma (-p)\, p(p-1)\cdots(p-n+1)2^n(2j_{n,p,a}-n-a)^{p-n}$$ for some real $j_{n,p,a}\in(0,n)$, so that $$i^{-p}(-1)^n I_{n,p,a}=\Gamma (n-p)\,2^n(2j_{n,p,a}-n-a)^{p-n}>0.$$ Quite similarly, $$i^p(-1)^n I_{n,p,a}=\Gamma (n-p)\,2^n(n+a-2j_{n,p,a})^{p-n}>0$$ if $a>n$ and $p$ is a non-integer in the interval $(-1,n)$. So, $I_{n,p,a}\ne0$ if $|a|>n$ and $p$ is a non-integer in the interval $(-1,n)$.
By using the l'Hospital rule as mentioned in the question statement, one should likely get the same result when $|a|>n$ and $p$ is an integer in the interval $(-1,n)$.
One can also consider the asymptotic behavior of $I_{n,p,a}$ as $n\to\infty$. Writing $e^{it}-e^{-it}$ (in the expression for $g_{n,a}(t)$) as $2i\sin t$ and then writing $\sin(\frac\pi2+k\pi+h)=(-1)^k\,e^{-h^2/(2+o(1))}$ for $k=0,1,\dots$ and $h\to0$, one sees that $$I_{n,p,a}\sim\frac{C_{p,a}}{\sqrt n}L(b/2,1/2,p+1)$$ as $n\to\infty$, where $C_{p,a}$ is manifestly nonzero and does not depend on $n$, $b=b_{n,a}$ equals $a$ or $a+1$ depending on whether $n$ is even or odd, and $$L(\lambda,c,s):=\sum_{k=0}^\infty\frac{\exp(2\pi i\lambda k)}{(k+c)^s}$$ is the Lerch zeta-function; see e.g. Wikipedia. It is now shown at positivity property of the Lerch zeta-function that $\Re L(b/2,c,p+1)>0$ for all real $b\ne0$, all real $c>0$, and all real $p>-1$. So, $I_{n,p,a}\ne0$ for large enough $n$. (If $b=0$, then $a\in\{-1,0\}$, in which case it is easy to see that $I_{n,p,a}\ne0$.)
