Divisibility property for sequence $a_{n+2}=-2(n-1)(n+3)a_n-(2n+3)a_{n+1}$

The sequence starts $$0, 1, -5, 25, -105, 105, 5355, \dots$$

We can observe that the statement is true not only for primes, but for odd numbers in general. Even though recurrences might be better for solving this kind of problem, here we can go for closed form formula (luckily there is one!). I have used approach inspired by http://mathforum.org/library/drmath/view/67314.html.

Currently the coefficients contain quadratic polynomials in $n$, but we can make them linear, which makes it a bit easier to work with. With substitution $b_n=\frac{a_n}{(n+2)!}$ and some algebra we get

$$(n+4)b_{n+2}+(2n+3)b_{n+1}+(2n-2)b_n=0$$

Now that the coefficients are linear, let's try to find its generating function $y(x)=\sum_{n \geq 1}b_n x^n$. By summing the whole equation and again some technical steps, we arrive at

$$ y'(x)(2x^3+2x^2+x)=\frac{1}{6}x^2+(2x^2-x-2)y(x) $$ which is a linear differential equation with variable coefficients, and can be solved to $$ y(x)=-\frac{1}{3x}-\frac{1}{9x^2}-\frac{5x}{18}-\frac{1}{2}+\frac{(2x^2+2x+1)^{3/2}}{9x^2} $$ (I have used CAS system for solving this differential equation, but it is just a technicality to show, also it is trivial to verify by differentiation).

Now to get information about the coefficients, let's expand it into the series. The power on the right is by Binomial series $$ (2x^2+2x+1)^{3/2}=\sum_{k=0}^{\infty}\binom{3/2}{k}2^k(x^2+x)^k $$ Using Binomial theorem for the inner sum and playing with the indicies, we obtain $$ [x^n](2x^2+2x+1)^{3/2} = \sum_{n/2 \leq k \leq n}\binom{3/2}{k}\binom{k}{2k-n}2^k $$ Dividing by $9x^2$ and noticing first few terms in result series are equal to $\frac{1}{9x^2}+\frac{1}{3x}+\frac{1}{2}+\frac{5x}{18}$, we can cancel those out in the $y(x)$, and so $$ b_n=\sum_{(n+2)/2 \leq k \leq n+2}\binom{3/2}{k}\binom{k}{2k-n-2}2^k $$ and in turn $$ a_n=(n+2)!\sum_{(n+2)/2 \leq k \leq n+2}\binom{3/2}{k}\binom{k}{2k-n-2}2^k $$ Looking at first few individual terms in a sum for couple values of $n$, we can observe those are all integers. This suggests we can further simplify the expression. Writing out the definitions and all the factorials, one eventually finds for $n>1$: $$ \boxed{a_n = \sum_{(n+2)/2 \leq k \leq n+2}\binom{k}{n+2-k}\frac{(n+2)!}{k!}\frac{(2k-5)!!}{3}(-1)^k} $$ We can see that for $k<n+2$ the terms above are divisible by $n+2$. So we have \begin{align} a_n &\equiv \binom{n+2}{0}\frac{(n+2)!}{(n+2)!}\frac{(2(n+2)-5)!!}{3}(-1)^{n+2}\\ &= \frac{(2n-1)!!}{3}(-1)^{n+2}\\ &= (2n-1)(2n-3)\cdots 5 \cdot (-1)^{n+2} \pmod{n+2} \end{align} If $n\geq 3$ is odd, the $n+2$ is odd as well, so the product above will contain it and thus will be divisible by it. In other words $a_n \equiv 0 \pmod {n+2}$, which is what we wanted to prove.

The claim that $a_{n+1} \equiv 0 \pmod {n+2}$ now follows for example by writing the original equation as

$$ a_{n+1}=-2(n-2)(n+2)a_{n-1}-(2n+1)a_n $$ which modulo $n+2$ yields $$ a_{n+1} \equiv 3a_n \equiv 0 \pmod {n+2} $$ I might have omitted few details here and there, but main idea is hopefully clear.

Define $$a_1=0,\; a_2=1,\;\text{and}\quad a_{n+2}=-2(n-1)(n+3)a_n-(2n+3)a_{n+1}.\tag{0}$$ Define $$\;R(x) := \sqrt{1-4x} = 1-2\sum_{n=1}^\infty C_{n-1}x^n = \sum_{n=0}^\infty {2n \choose n}\frac{x^n}{1-2n}\;$$ which is the generating function of the OEIS integer sequence A002420. Define $$\;B(x) := \frac{ (1+4x+8x^2)R(-x-2x^2)-(1+6x+18x^2+20x^3)}{2x^2}$$ $$ = 3x^2 - 6x^3 + 10x^4 - 12x^5 + 3x^6 + \dots =: \sum_{n=0}^\infty b_nx^n$$ which is the generating function for $b_n.\;$ It is easy to check that $b_n$ is also an integer sequence and $\;b_1=0, b_2=3.\;$ Using derivatives we get the differential equation $$\;(8x^2-2x-2)B(x)+12x^2=x(1+4x+8x^2)\frac{d}{dx}B(x)\;$$ and using some algebra on power series coefficients we get the equation $$\;(2+n)b_{n}=(2-4n)b_{n-1}+(24-8n)b_{n-2}\quad\forall\;n>2.$$

Redefine $$\;a_n:=\frac{(n+2)!}{2^{n+1} 9}b_n.\tag{1}$$ It is easy to check that this satisfies the same initial values and recursion as given in $(0)$, but is it an integer sequence? It will be if $\;2^{n+1}\;|\;(n+2)!\;b_n\;$ for all $n>1.$ The key results to prove is that if $\;C_n\;$ is the Catalan numbers and if $\;v_2(n)\;$ is the 2-adic valuation of $n$, then

$$v_2(b_n) = v_2(C_{n+1}) = v_2\Big(\frac{2^{n+1}}{(n+2)!}\Big),\; v_2(a_n)=0\quad \forall n>1. \tag{2}$$

Now the original equation $(0)$ already implies that $\;a_n\;$ is an integer sequence using induction, $\;a_n\;$ is odd if $n>1$, and the rest of equation $(2)$ is true.

However, there is an alternative approach using the exponential generating function of $\;a_n,\;$ namely $\;A(x):=\sum_{n=0}^\infty a_nx^n/n!.\;$ Using equation $(0)$ for $\;a_n\;$ we get the differential equation $$ 1 = -6A(x) +(3+6x)\frac{d}{dx}A(x) + (1+2x+2x^2)\frac{d^2}{dx^2}A(x).$$ Taking the derivative of this equation gives another differential equation $$ 0 = 5(1+2x)\frac{d^2}{dx^2}A(x) +(1+2x+2x^2)\frac{d^3}{dx^3}A(x).$$ Now define the sequence $\;y_n(x)\;$ by $\,y_0(x) := 1+O(x),\,$ and recursion $$ y_{n+1}(x) := 1-\int_0^x 5(1+2x)/(1+2x+2x^2)\;y_n(x)\; dx. \tag{3}$$ The limit as $\;n\to\infty$ of $\;y_n(x) =\frac{d^2}{dx^2}A(x).\;$ It is easy to check that the recursion $(3)$ preserves exponential generating functions of integer sequences. Thus, again $\;a_n\;$ is an integer sequence.

If $n>0$ is odd, then $(n+2)\;|\;a_n$ and if $n>2$ is even, then $(n+1)\;|\;a_n.\;$ This since one of the $(n\!+\!1)(n\!+\!2)$ factors in the numerator of $(1)$ are not cancelled by the $2^{n+1}$ in the denominator.