For $q \in \mathbb{N}$, consider the following infinite product:
$$\Pi_q(x) = \prod_{m=1}^\infty \left(1-\frac{x^q}{n^q}\right) = \prod_{r=1}^q \frac{1}{\Gamma(1-\zeta_q^r x)}, \qquad (1)$$
Last equality follows from Weierstrass factorization theorem, together with the observation that:
$$1-\frac{x^q}{n^q} = \left( 1 - \frac{\zeta_q x}{n}\right)\left( 1 - \frac{\zeta_q^2 x}{n}\right)\dots\left( 1 - \frac{\zeta_q^q x}{n}\right) = \prod_{r=1}^q \left( 1 - \frac{\zeta_q^r x}{n}\right) $$
with $\zeta_q$ the $q$th root of unity given by $\zeta_q= \exp(2\pi i/q)$ . For even $q$, it is possible to express the product of gamma functions as a closed formula, thanks to Euler's reflection formula:
$$\Gamma(1-z) \Gamma(1+z) = \frac{\pi z}{\sin \pi z}, \qquad z\in \mathbb{C}-\mathbb{Z}$$
Then for $q$ even, $\Pi_q(x)$ in $(1)$ can be written as:
$$\Pi_q(x) = \prod_{r=1}^q \frac{1}{\Gamma(1-\zeta_q^r x)} = \prod_{r=1}^{q/2} \frac{1}{\Gamma(1-\zeta_q^r x)\Gamma(1+\zeta_q^r x)} = \prod_{r=1}^{q/2} \frac{\sin(\pi\zeta_q^r x)}{\pi \zeta_q^{r} x}, \qquad (2)$$
because, for $q$ even, if $\zeta_q$ is a root of unity, also $-\zeta_q$ is a root of unity. Some specific cases of the formula are:
$$ \Pi_2(x) = \frac{\sin\pi x}{\pi x} = 1 - \underbrace{(\pi^2/6)}_{\zeta(2)}x^2 + \dots $$
$$ \Pi_4(x) = \frac{\sin\pi x \sinh \pi x}{\pi^2 x^2} = 1-\underbrace{(\pi^4/90)}_{\zeta(4)}x^4 + \dots $$ $$ \Pi_6(x) = \frac{\sin\pi x(\cos\pi x - \sqrt{3}\cosh\pi x) }{\pi^3 x^3} $$ $$ \Pi_8(x) = \frac{\sin\pi x \sinh \pi x\left[\cos^2(\pi x/\sqrt{2})\sinh^2(\pi x/\sqrt{2}) + \sin^2(\pi x/\sqrt{2})\cosh^2(\pi x/\sqrt{2})\right] }{\pi^4 x^4}$$
My question has to do with the general formula for $\Pi_q(x)$ when $q$ is even in $(2)$. From the above instances of $\Pi_q(x)$, it seems that the numerator of the formulas that define $\Pi_q(x)$ looks like a function with cyclic derivatives (aside one factor). That is, if we define the function $\sigma^n_0(x)$ as the solution of the differential equation: $$\frac{\text{d}^n\sigma^n_0(x)}{\text{d}x^n} + \sigma^n_0(x) = 0 $$
with $\sigma^n_0 (0) =1, D^k\sigma^n_0(0) =0$ for $1\le k \le n-1$. Then, we define in addition: $$\sigma^n_k(x) = \frac{\text{d}^k\sigma^n_0(x)}{\text{d}x^k}, \text{ for } 1\le k \le n-1$$ Some examples of $\sigma^{2q}_q(x)$ are: $$\sigma^2_1(x) = \sin(x)$$ $$\sigma^4_2(x) = -\sin\left(\frac{x}{\sqrt(2)}\right) \sinh\left(\frac{x}{\sqrt(2)}\right)$$ It is clear that for $q=2,4,6$ and $8$, the numerators in the formulas for $\Pi_q(x)$ contain $\sigma^{2q}_q(x)$. But we need to show that the $q$ derivative of:
$$\frac{\text{d}^q}{\text{d}x^1} \left( \prod_{r=1}^{q/2} \sin(\pi\zeta_q^r x) \right) \propto \prod_{r=1}^{q/2} \sin(\pi\zeta_q^r x)$$
that is, the $n$th derivative of the function is proportional to the function itself. My attempt was to use multinomial theorem for derviatives:
$$D^{2m} (\Phi_1(x)\dots\Phi_m(x)) = \sum_{k_1+k_2+\cdots+k_m=2m} {2m \choose k_1, k_2, \ldots, k_m} D^{k_1}\Phi_1 D^{k_2}\Phi_2 \dots D^{k_m}\Phi_m $$
And then use that an even-ordered derivative ($2p$) of $\sin(\zeta_q^r x)$ is proportional to $(-1)^p(\zeta_q^{2p}) \sin(\zeta_q^r x)$, while an odd-order derivative ($2p+1$) is proportional to $(-1)^{2p} (\zeta_q^{2p+1}) \cos(\zeta_q^r x)$. And then, looking at the structure of the multinomial coefficients, check the cancellations of terms with $\cos$. Although the approach is very elementary, the task of checking the cancellations well seems cumbersome and tiring. Is there a simpler way to prove this?
Note that once a closed form for $\Pi_q(x)$ has been determined, the sum of reciprocal powers can be calculated as:
$$\zeta(k) = \frac{1}{1^k} + \frac{1}{2^k} + \frac{1}{3^k} + \dots = \frac{1}{k!} \frac{\text{d}^k \Pi_k(0)}{\text{d}x^k}$$
Although the problem is to find a closed-form expression for $\Pi_q(x)$ in terms of elementary functions (if it exists). For $q$ even, Euler proved around 1735 this was always possible. Unfortunately, for $p = 3$ or $p=2k+1$, the function $\Pi_3(x)$ does not seem expressible in terms of elementary functions, so no analytical method of determining the Apéry's constant $\zeta(3)$ is known: $$ \Pi_3(x) = \frac{1}{\Gamma(1-x)\Gamma(1-\frac{-1+i\sqrt{3}}{2}x)\Gamma(1-\frac{-1-i\sqrt{3}}{2}x)} \approx 1 -\zeta(3)x^3 + \frac{\zeta(3)^2-\zeta(6)}{2} x^6 + \dots $$