It can be shown that $$\int^\pi_{-\pi}(\text{coc }x)^a dx=\frac{2\pi\Gamma(a+1)}{\Gamma^2\left(\frac a2+1\right)}$$
Differentiation under the integral sign and the identity
$$\psi_n(1)=(-1)^{n+1}n!\zeta(n+1)$$ together lead to a clear relationship between $c_n$ and $\zeta$ function.
I will elaborate very soon.
Real analysis derivation:
Recall that $$\mathcal B(x,y)=2\int^{\pi/2}_{0}\sin^{2x-1}t \cos^{2y-1}t dt$$
Therefore,
$$\begin{align}
\int^\pi_{-\pi}(\text{coc }x)^a dx
&=2^a\cdot 2\int^\pi_{0}\cos^a\left(\frac x2\right) dx \\
&=2^{a+1}\cdot 2\int^{\pi/2}_0\cos^a (u) du \\
&=2^{a+1}\mathcal B\left(\frac12,\frac{a+1}2\right) \\
&=2^{a+1}\frac{\sqrt \pi }{\Gamma\left(\frac a2 +1\right)}\cdot \Gamma\left(\frac{a+1}2\right) \\
&=2^{a+1}\frac{\sqrt \pi }{\Gamma\left(\frac a2 +1\right)}\cdot \frac{2^{1-(a+1)}\sqrt\pi\cdot\Gamma(a+1)}{\Gamma\left(\frac{(a+1)+1}{2}\right)} \qquad{(\star)}\\
&=\frac{2\pi\Gamma(a+1)}{\Gamma^2\left(\frac a2+1\right)}
\end{align}
$$
$(\star):$ Legendre duplication formula $$\Gamma\left(\frac v2\right)=\frac{2^{1-v}\sqrt\pi\cdot\Gamma(v)}{\Gamma\left(\frac{v+1}{2}\right)}$$ is used.
Complex analysis derivation:
$$\begin{align}
\int^{\pi}_{-\pi}(\text{coc }x)^a dx
&=\int^{\pi}_{-\pi}(e^{ix/2}+e^{-ix/2})^a dx \\
&\stackrel{z=e^{ix}}{=}\oint_{|z|=1}\left(\sqrt z+\frac1{\sqrt z}\right)^a\frac{dz}{iz} \\
&=\oint_{|z|=1}\frac{(z+1)^a}{iz^{a/2+1}}dz \\
&=-\lim_{\epsilon\to0^+}\left(\int^{0+i\epsilon}_{-1+i\epsilon}+\int^{-1-i\epsilon}_{0-i\epsilon}\right)\frac{(z+1)^a}{iz^{a/2+1}}dz \\
&=i\left(\int^0_{-1} \frac{(z+1)^a}{e^{i\pi(a/2+1)}|z|^{a/2+1}}dz+\int^{-1}_0 \frac{(z+1)^a}{e^{-i\pi(a/2+1)}|z|^{a/2+1}}dz \right) \qquad(1)\\
&=i\left(-e^{-i\pi a/2}\int^1_{0} \frac{(1-z)^a}{z^{a/2+1}}dz+e^{i\pi a/2}\int^{1}_0 \frac{(1-z)^a}{z^{a/2+1}}dz \right) \\
&=-2\sin\frac{\pi a}{2}\mathcal B\left(a+1,-\frac a2\right) \\
&=\frac{2\pi\Gamma(a+1)}{\Gamma^2\left(\frac a2+1\right)} \qquad{(2)}\\
\end{align}
$$
$(1)$: Consider a key hole contour avoiding the principal logarithmic branch cut on the negative real axis, and apply Cauchy's integral theorem (no singularities are enclosed).
$(2)$: By applying the Gamma reflection formula.
Hence,
$$c_n = \left(\frac{\partial}{\partial a}\right)^n \underbrace{\frac{2\Gamma(a+1)}{\Gamma^2\left(\frac a2+1\right)}}_{f(a)}\bigg\vert_{a=0} \\$$
For example, differentiating twice gives
$$\begin{align}
c_2
&=\frac{2\Gamma(a+1)}{\Gamma^2\left(\frac a2+1\right)}\left[\psi_0(a+1)-\psi_0\left(1+\frac a2\right)\right]^2+\frac{2\Gamma(a+1)}{\Gamma^2\left(\frac a2+1\right)}\left(\psi_1(a+1)-\frac12\psi_1\left(1+\frac a2\right)\right)\bigg\vert_{a=0} \\
&=2\cdot\frac12\psi_1(1) \\
&=(-1)^{1+1}1!\zeta(1+1) \\
&=\zeta (2)
\end{align}
$$
Other $c_n$ can be found similarly. I tried calculating $c_5$ by hand, and it turns out the algebra is pretty tedious.
It it useful to define $p_k=\psi_k(a+1)-\frac1{2^k}\psi_k\left(\frac a2+1\right)$, since
$$f’=fp_0\qquad\qquad p_n’=p_{n+1}$$
Doing some algebra, I got
$$\frac{f^{(5)}}{f}=p_0^5+10p_0^3p_1+15p_0p_1^2+10p_0^2p_2+10p_1p_2+5p_0p_3+p_4$$
As $p_0(0)=0$,
$$c_5=2(10p_1p_2+p_4)=20\cdot\frac{\pi^2}{12}\frac{-3\zeta(3)}2+2\cdot\frac{-45\zeta(5)}2$$
$$\implies c_5=-\frac52\pi^2\zeta(3)-45\zeta(5)$$
It can be seen that the numerical factors grow quite quickly. One may be interested in the asymptotics of $c_n$.
Notice that $\frac{c_n}{n!}$ is the $n$th coefficient of the Maclaurin series of $f$. Due to the nearest pole at $a=-1$, $$\frac{|c_{n+1}/(n+1)!|}{|c_n/n!|}\sim 1\implies |c_{n+1}|\sim (n+1)|c_n|$$ implying factorial growth.
$f(a)$ satisfies
$$f'(a)=f(a)\underbrace{\left((\psi_0(a+1)-\psi_0\left(\frac a2+1\right)\right)}_{\gamma_1(a)}$$
In general $f^{(n)}(a)=f(a)\gamma_n(a)$ where $$\gamma_{n+1}=\gamma_1\gamma_n+\gamma_{n}'$$
Since $\gamma_1(0)=0$,
$$c_n=2\gamma_n'(0)$$
