Define the triple factorial, $n!!!$, as a continuous function for $n \in \mathbb{C}$

Question

For anyone unfamiliar with multifactorial notation, I will give a quick rundown of it (at least, to the best of my understanding) for non-negative integer values of $n$:

$$n!=n(n-1)(n-2)(n-3)...(n-a), (n-a) > 0$$ $$n!!=n(n-2)(n-4)(n-6)...(n-a), 2 \geq (n-a) > 0$$ $$n!!!=n(n-3)(n-6)(n-9)...(n-a), 3 \geq (n-a) > 0$$

A more generalised expression can be given like so, where $k$ represents the number of factorial symbols:

$$n!^{k}=\left( \begin{cases} 1 & n=0\\ n & 1\leq n\leq k \\ n(n-k)(n-2k)(n-3k)...(n-a) & n>k \end{cases} \right), k \geq (n-a) > 0$$

Now, this is great for when you're working with (primarily) positive integers, but I'm curious how you'd go about extending the definition such that it will be valid for all real and complex numbers. While a general definition for any multifactorial would be amazing, I am primarily just looking for for the definition regarding the triple factorial.

I already know it's possible to do so for the double factorial; $z!!=2^{(1+2z-\cos(\pi z))/4}\pi^{(\cos(\pi z)-1)/4}\Gamma(z/2+1), z \in \mathbb{C}$, so it has me hopeful that it's also possible to define $z!!!$ in a similar such manner.

Answer 1

The expression mentioned in your question is a reflection of the fact that there are actually two "nice" extensions of the double factorial (one agreeing with it on even numbers and the other agreeing on odd numbers) and you can only combine them in a somewhat ad-hoc way.

Namely, the even and odd double factorials are

$$z!!_0 = \frac{2^{\frac{z}{2}} \Gamma\left(1+\frac{z}{2}\right)}{\Gamma\left(1+\frac{0}{2}\right)} = 2^{\frac{z}{2}} \Gamma\left(1+\frac{z}{2}\right)$$

and

$$z!!_1 = 1 \frac{2^{\frac{z-1}{2}} \Gamma\left(1+\frac{z}{2}\right)}{\Gamma\left(1+\frac{1}{2}\right)} = \sqrt{\frac{2}{\pi}} 2^{\frac{z}{2}} \Gamma\left(1+\frac{z}{2}\right).$$

The expression in your question is equal to

$$z!! = T_2(z) \: 2^{\frac{z}{2}} \Gamma\left(1+\frac{z}{2}\right),$$

where $T_2(z) = (\sqrt{2/\pi})^{(1-\cos (\pi z))/2}$ is equal to $1$ when $z$ is an even integer and to $\sqrt{2/\pi}$ when $z$ is an odd integer. I say that this is an ad-hoc construction because there are clearly many more possible $T_2(z)$ with this property (for example $(\sqrt{2/\pi})^{(1-\cos (13\pi z))/2}$).

Similarly, for the triple factorial we have three "nice" extensions, based on the triple factorials of numbers congruent to $0, 1, 2 \pmod 3$:

$$z!!!_0 = \frac{3^{\frac{z}{3}} \Gamma\left(1+\frac{z}{3}\right)}{\Gamma\left(1+\frac{0}{3}\right)},$$

$$z!!!_1 = \frac{3^{\frac{z-1}{3}} \Gamma\left(1+\frac{z}{3}\right)}{\Gamma\left(1+\frac{1}{3}\right)},$$

$$z!!!_2 = 2 \frac{3^{\frac{z-2}{3}} \Gamma\left(1+\frac{z}{3}\right)}{\Gamma\left(1+\frac{2}{3}\right)}.$$

So an extension of the triple factorial to $\mathbb{C}$ would be

$$z!!! = T_3(z) \: 3^{\frac{z}{3}} \Gamma\left(1+\frac{z}{3}\right),$$

where $T_3(z)$ is any function interpolating the three constants $1$, $3^{\frac{-1}{3}}/\Gamma\left(1+\frac{1}{3}\right)$ and $2 \cdot 3^{\frac{-2}{3}}/\Gamma\left(1+\frac{2}{3}\right)$ for $\mathbb{Z}\ni z=0,1,2 \pmod 3$ respectively, for example

$$T_3(z) = \left(\frac{3^{\frac{-1}{3}}}{\Gamma\left(1+\frac{1}{3}\right)} \right) ^{(1+2\cos(2\pi (z-1)/3))/3} \left(2 \frac{3^{\frac{-2}{3}}}{\Gamma\left(1+\frac{2}{3}\right)} \right)^{(1+2\cos(2\pi (z-2)/3))/3}.$$

You can see the resulting function here.

---

We can argue similarly for the higher-order multifactorials $z!^k$. The final expression is

$$z!^k = T_k(z) \: k^{\frac{z}{k}} \Gamma\left(1+\frac{z}{k}\right),$$

and a possible $T_k(z)$ generalizing the expressions for $T_2(z)$ and $T_3(z)$ given above is

$$T_k(z) = \prod_{j=1}^{k} \left(j\frac{k^{-j/k}}{\Gamma\left(1+\frac{j}{k}\right)} \right)^{E(k,j;z)},$$

where the exponent $E(k,j;z)=\frac{1}{k}\sum_{l=1}^k \cos\left(2\pi l \frac{(z-j)}{k}\right)$ is $1$ when $z = j \pmod k$ and vanishes when $z$ is any other integer. You can see its graph here (the slider also seems to allow fractional values of $k$, but I wouldn't trust that case).

Answer 2

EDIT 29.12.19

After a more thorough revision I conclude that text below is not to be considered as a solution.

Please see paragraph 3. in the discussion where my doubts started.

Here is a comparision of my approach (W) changed to standard initial conditions and that of pregunton (P). I use $z$ as the argument hencforth.

1. both W and P satisfy the recursion $(1)$ even for non integer values of $z$
2. P is real for real $z$, W is complex for non integer real $z$
3. P and Re(W) differ for non integer real $z$

Hence I conclude that P has given a satisfactory solution, and that W should be dismissed.

The trouble arose because I had erroneously thought the solution of the recursion to be unique.

Maybe my approach can be saved requesting $a_k(z)$ be real on the real $z$-axis.

Abstract

This solution for the $k$-fold faculty is in principle complete and straightforward. It combines linear recursion, roots of unity and linear algebra. Examples for $k=1..5$ are given. Hints to the OEIS database are given for the (sub-)sequences encountered.

Derivation of the solution

Let $a_{k}(n)$, the $k$-fold faculty of $n$ be defined by the recursion

$$a_{k} (n) = n a_{k} (n-k)\tag{1}$$

with the inital conditions

$$a_{k}(m) = 1, m=0,1,...,k-1\tag{2}$$

Now, since $(1)$ is a linear homogeneous difference equation of order $k$ its solution is a linear combination of $k$ linearly independent solutions.

We can find a particular solution letting with some parameter $\lambda$

$$a_{k} (n) = \lambda^\frac{n}{k} b_{k}(1+\frac{n}{k})\tag{3}$$

Inserting this in $(1)$ we have

$$ \lambda^\frac{n}{k} b_{k}(1+\frac{n}{k})= n \lambda^\frac{n-k}{k} b_{k}(1+\frac{n-k}{k})$$

or

$$ b_{k}(1+\frac{n}{k})= \frac{n}{\lambda} b_{k}(\frac{n}{k})\tag{4}$$

Letting $\lambda=k$ the simplified recursion $(4)$ converts into the defining recursion for the $\Gamma$-function and we have

$$a_{k} (n)=A k^{\frac{n}{k}} \Gamma(1+\frac{n}{k})\tag{5}$$

with some constant $A$.

In order to find more solutions we employ the method of varying the "constant" writing

$$a_{k} (n)=A(n)k^{\frac{n}{k}}\Gamma(1+\frac{n}{k})\tag{6}$$

Inserting this into $(5)$ gives the recursion

$$A(n) = A(n-k)\tag{7a} $$

The $k$ linearly independent solutions of which are roots of unity

$$A(n,m) = e^{2 \pi i n \frac{m}{k}}, m=0,1,...,(k-1)\tag{7b}$$

Hence, with constants $c_m$, we can write down the general solution of $(1)$ as

$$a_{k}(n) = k^{\frac{n}{k}} \Gamma(1+\frac{n}{k}) \sum_{m=0}^{k-1} c_{m} e^{2 \pi i m \frac{n}{k}}\tag{8}$$

The $k$ constants $c_{m}$ have to be determined by the $k$ initial conditions $(2)$.

Examples

While $(8)$ together with $(2)$ give the complete solution to the problem it is still more or less implicit, and therefore we calculate explit formulae for some small values of $k$.

$k=1$

We have $a_1(n)=c_0 \Gamma(1 + n)$ and from $(2)$ follows $c_0=1$ so that

$$a_1(n) = \Gamma(1+n) = n!\tag{e1}$$

as expected.

$k=2$

We have $a_2(n) = 2^{n/2} \left(c_0+c_1 e^{i \pi n}\right) \Gamma \left(\frac{n}{2}+1\right)$.

The initial conditions lead to the equations $1=c_0+c_1, 1 = \sqrt{2} (c_0 - c_1)$ giving $c_0= \frac{1}{2}\left(1+\sqrt{\frac{2}{\pi}}\right),c_1 =\frac{1}{2}\left(1-\sqrt{\frac{2}{\pi}}\right)$.

Finally

Real $n$:

$$a_2(n) = 2^{\frac{n}{2}}\Gamma \left(\frac{n}{2}+1\right) \frac{1}{2}\left(\left(1+\sqrt{\frac{2}{\pi }}\right)+\left(1-\sqrt{\frac{2}{\pi }}\right) e^{i \pi n}\right)\tag{e2a} $$

Graph:

Integer $n$:

For integer we have the simplification $n$ $e^{i \pi n}=(-1)^n$.

Since for even $n$: $\frac{1}{2}\left(\left(1+\sqrt{\frac{2}{\pi }}\right)+\left(1-\sqrt{\frac{2}{\pi }}\right) (-1)^n\right) = 1$ for even $n$ and $=\sqrt{\frac{2}{\pi }}$ for odd $n$ we can also write

$$a_2(n) = 2^{\frac{n}{2}}\Gamma \left(\frac{n}{2}+1\right)(\sqrt{\frac{2}{\pi }})^{\frac{1-(-1)^n}{2}}\tag{e2b} $$

A similar form was given in the OP.

The first few values are

$$\{a_2(n) \}_{n=0}^{n=11} = \{1,1,2,3,8,15,48,105,384,945,3840\}$$

This sequence is listed here https://oeis.org/A006882 as double factorials n!!: a(n) = n*a(n-2) for n > 1, a(0) = a(1) = 1.

$k=3$

This case is even more interesting since $n$ has three values $\mod(3)$

We have, to begin with,

$$a_{3}(n) = 3^{n/3} \left(c_0+c_1 e^{\frac{2 i \pi n}{3}}+c_2 e^{\frac{4 i \pi n}{3}}\right) \Gamma \left(\frac{n}{3}+1\right)\tag{e3a}$$

The inital conditions lead to these linear equations for the coefficients

$$ \begin{array}{c} 1=c_0+c_1+c_2 \\ 1=\sqrt[3]{3} \left(c_0+e^{\frac{2 i \pi }{3}} c_1+e^{-\frac{1}{3} (2 i \pi )} c_2\right) \Gamma \left(\frac{4}{3}\right) \\ 1=3^{2/3} \left(c_0+e^{-\frac{1}{3} (2 i \pi )} c_1+e^{\frac{2 i \pi }{3}} c_2\right) \Gamma \left(\frac{5}{3}\right) \\ \end{array} $$

Real $n$:

The solution of which gives rise to

$$a_3(n) = 3^{n/3} \Gamma \left(\frac{n}{3}+1\right)\left( \frac{1}{3}\left(1+2 (-1)^n \cos \left(\frac{\pi n}{3}\right)\right) \\ +\frac{1}{3} \left(1+(-1)^{2/3} e^{\frac{2 i \pi n}{3}}-\sqrt[3]{-1} e^{\frac{4 i \pi n}{3}}\right) \left(\frac{3^{5/6} \Gamma \left(\frac{1}{3}\right)}{ 4 \pi }\right) \\ +\frac{1}{3}\left((-1)^{n+1} \left(\sqrt{3} \sin \left(\frac{\pi n}{3}\right)+\cos \left(\frac{\pi n}{3}\right)\right)+1\right)\left(\frac{6 \sqrt[6]{3} \Gamma \left(\frac{2}{3}\right)}{ 4 \pi }\right)\right)\tag{e3b} $$

Graph:

Integer $n$:

Simplifying according to the remainder$\mod(3)$ we have

$a_{3}(n) = \left\{ \begin{array}{ll} 3^j \Gamma (j+1) & n=3j + 0 \; \;\;\{\}_{j=0}^{j=5}= \{1,3,18,162,1944,29160\},\; \text{https://oeis.org/A032031}\\ \frac{3^{j+\frac{3}{2}} \Gamma \left(\frac{2}{3}\right) \Gamma \left(j+\frac{4}{3}\right)}{2 \pi } & n=3j+1 \;\;\; \{\}_{j=0}^{j=5} = \{1,4,28,280,3640,58240\},\; \text{https://oeis.org/A007559}\\ \frac{3^{j+\frac{3}{2}} \Gamma \left(\frac{1}{3}\right) \Gamma \left(j+\frac{5}{3}\right)}{4 \pi } & n=3j+2 \;\; \; \{\}_{j=0}^{j=5}=\{1,5,40,440,6160,104720\}, \text{https://oeis.org/A034000}\\ \end{array} \right.\tag{e3c} $

$k=4$

Real $n$:

$$a_{4}(n) = \frac{2^{\frac{n}{2}-2}}{3\pi} \Gamma \left(\frac{n}{4}+1\right) \left(-3 \sqrt{\pi } \left(-1+e^{\frac{i \pi n}{2}}\right) \left(1+e^{i \pi n}\right)+3 \pi \left(1+e^{\frac{i \pi n}{2}}\right) \left(1+e^{i \pi n}\right)-\left(-1+e^{i \pi n}\right) \left(4 i e^{\frac{i \pi n}{2}} \left(\Gamma \left(\frac{5}{4}\right)-2 \Gamma \left(\frac{7}{4}\right)\right)+\Gamma \left(\frac{1}{4}\right)+8 \Gamma \left(\frac{7}{4}\right)\right)\right)$$

Graph:

Integer $n$:

$a_{4}(n) = \left\{ \begin{array}{ll} 4^j \Gamma (j+1) & n=4j + 0 \; \;\;\{\}_{j=0}^{j=5}= \{1,4,32,384,6144,122880\},\; \text{https://oeis.org/A047053}\\ \frac{2^{2 j+\frac{7}{2}} \Gamma \left(\frac{7}{4}\right) \Gamma \left(j+\frac{5}{4}\right)}{3 \pi } & n=4j+1 \;\;\; \{\}_{j=0}^{j=5} = \{1,5,45,585,9945,208845\},\; \text{https://oeis.org/A007696}\\ \frac{2^{2 j+1} \Gamma \left(j+\frac{3}{2}\right)}{\sqrt{\pi }} & n=4j+2 \;\; \; \{\}_{j=0}^{j=5}=\{1,6,60,840,15120,332640\}, \text{https://oeis.org/A000407}\\ \frac{2^{2 j+\frac{3}{2}} \Gamma \left(\frac{1}{4}\right) \Gamma \left(j+\frac{7}{4}\right)}{3 \pi } & n=4j+3 \;\; \; \{\}_{j=0}^{j=5}=\{1, 7, 77, 1155, 21945, 504735\}, \text{https://oeis.org/A034176}\\ \end{array} \right. $

$k=5$

Real $n$: lengthy expression, dropped here

Graph:

Integer $n$:

$a_{5}(n) = \left\{ \begin{array}{ll} 5^j \Gamma (j+1) & n=5j + 0 \; \;\;\{\}_{j=0}^{j=5}= \{1,5,50,750,15000,375000\},\; \text{https://oeis.org/A052562} \\ \frac{\sqrt[10]{\frac{1}{2} \left(25-11 \sqrt{5}\right)} \;5^{j+\frac{11}{5}} \Gamma \left(\frac{9}{5}\right) \Gamma \left(j+\frac{6}{5}\right)}{8 \pi } & n=5j+1 \;\;\; \{\}_{j=0}^{j=5} = \{1,6,66,1056,22176,576576\},\; \text{https://oeis.org/A008548} \\ \frac{\sqrt[10]{\frac{1}{2} \left(25+11 \sqrt{5}\right)}\; 5^{j+\frac{11}{5}} \Gamma \left(\frac{8}{5}\right) \Gamma \left(j+\frac{7}{5}\right)}{12 \pi } & n=5j+2 \;\; \; \{\}_{j=0}^{j=5}=\{1,7,84,1428,31416,848232\}, \text{https://oeis.org/A034323} \\ \frac{5^j \Gamma \left(j+\frac{8}{5}\right)}{\Gamma \left(\frac{8}{5}\right)} & n=5j+3 \;\; \; \{\}_{j=0}^{j=5}=\{1,8,104,1872,43056,1205568\}, \text{https://oeis.org/A034300} \\ \frac{\sqrt[10]{\frac{1}{2} \left(25-11 \sqrt{5}\right)}\; 5^{j+\frac{16}{5}} \Gamma \left(\frac{11}{5}\right) \Gamma \left(j+\frac{9}{5}\right)}{48 \pi } & n=5j+4 \;\; \; \{\}_{j=0}^{j=5}=\{1,9,126,2394,57456,1666224\}, \text{https://oeis.org/A034301} \\ \end{array} \right. $

Discussion

1. Simpler expressions for some $a_{k}(n)$ in terms of the Pochhammer symbol have been given in OEIS by G. C. Greubel, Aug 23 2019.

2. The expressions given here lead somewhat miraculously to integer numbers. The miracle can be solved by using the basic formula for the Gamma function $z \Gamma(z) = \Gamma(z+1)$ and the relation $\Gamma (z) \Gamma (1-z)=\pi \csc (\pi z)$.

3. I compared the results of the two solutions given up to now. For $k=3$, I found that they differ. That of pregunton (https://math.stackexchange.com/a/3488935/198592) agrees with OEIS A007661 "Triple factorial numbers a(n) = n!!!, defined by a(n) = n*a(n-3), a(0) = a(1) = 1, a(2) = 2. Sometimes written n!3." I call this the standard solution. Mine is different and not listed in OEIS. The reason is the different initial conditions. I had assumed that $a_3(2) = 1$, the "standard" is $a_3(2)=2$. For $k=2$ the standard is $a_2(0)=a_2(1)=1$, and consequently my series coincides with the standard. I don't quite understand at the moment why the generalization of the cases $k=1$ and $k=2$ should be a deviation from the constant initial conditions for all values below $k$ (my conditions $(2)$). And what which conditions would be "natural" for $k=4$ and higher $k$? I have found the answer here http://mathworld.wolfram.com/Multifactorial.html. The rule for the initial conditions is simply $a_{m}= m$ for $m=1..k+1$.

Answer 3

Solution as of 31.12.19

Summarizing statement

The two solutions proposed here (P by pregunton, and W by myself) follow two different strategies and lead to different solutions

P: calculates the series elements for each $mod(n,k)=0,1,...,(k-1)$ separately and then uses a switching function which selects each of these series from a general expression. The continuation to real (and complex) values of $n$ is then performed for the switching function.

W: here the basic starting point is the solution of the recursion relation which matches the initial conditions. The solution is tken to be a finite Fourier sum.

The results of the two approaches do not concide apart from integer arguments.

My reworked solution

The main differences are

(1) the intial conditions are now adapted to the standard values

(2) the solution to $$A(n) = A(n-k)$$ is now based on a real Fourier sum. It leaves some free parameters at our discretion which can be used to adapt the solution to a more specific one as e.g. to that of pregunton.

(3) the solution is now real on the real axis

Abstract

This solution for the $k$-fold faculty is in principle complete and straightforward. It combines linear recursion, Fourier sums and linear algebra. Examples for $k=1..5$ are given. Hints to the OEIS database are given for the (sub-)sequences encountered.

Derivation of the solution

Let $a_{k}(n)$, the $k$-fold factorial of $n$ be defined by the recursion

$$a_{k} (n) = n a_{k} (n-k)\tag{1}$$

with the inital conditions

$$a_{k}(0) = 1, a_{k}(m) = m, m=1,...,k-1\tag{2}$$

Now, since $(1)$ for integer $n$ is a linear homogeneous difference equation of order $k$ its solution is a linear combination of $k$ linearly independent solutions.

We can find a particular solution letting with some parameter $\lambda$

$$a_{k} (n) = \lambda^\frac{n}{k} b_{k}(1+\frac{n}{k})\tag{3}$$

Inserting this in $(1)$ we have

$$ \lambda^\frac{n}{k} b_{k}(1+\frac{n}{k})= n \lambda^\frac{n-k}{k} b_{k}(1+\frac{n-k}{k})$$

or

$$ b_{k}(1+\frac{n}{k})= \frac{n}{\lambda} b_{k}(\frac{n}{k})\tag{4}$$

Letting $\lambda=k$ the simplified recursion $(4)$ converts into the defining recursion for the $\Gamma$-function and we have

$$a_{k} (n)=A k^{\frac{n}{k}} \Gamma(1+\frac{n}{k})\tag{5}$$

with some constant $A$.

In order to find more solutions we employ the method of varying the "constant" writing

$$a_{k} (n)=A(n)k^{\frac{n}{k}}\Gamma(1+\frac{n}{k})\tag{6}$$

Inserting this into $(5)$ gives the recursion

$$A(n) = A(n-k)\tag{7a} $$

Obviously any function with period $k$ solves $(7a)$.

When it comes to periodic functions the most natural choice would be a Fourier series. We restrict the choice here to a finite Fourier sum:

$$A(n,k) =\sum_{m=0}^{k} c_{m} \cos(2 \pi m \frac{n}{k}) +s_0+\sum_{m=1}^{k} s_{m} \sin(2 \pi m \frac{n}{k})\tag{7b}$$

This expression contains $2k$ parameters. Since $k$ of them are defined by the initial conditions, $k$ are left undefined, and can be chosen arbitraryly. This reflects the fact that the extrapolation of a function defined for positive integers to real and complex values is to a high degree arbitrary.

Hence we can write down the general solution of $(1)$ as

$$a_{k}(n) = k^{\frac{n}{k}} \Gamma(1+\frac{n}{k})A(n,k)\tag{8}$$

(to be continued)

Examples

While $(8)$ together with $(2)$ give the complete solution to the problem it is still more or less implicit, and therefore we calculate explit formulae for some small values of $k$.

$k=1$

We have $a_1(n)=c_0 \Gamma(1 + n)$ and from $(2)$ follows $c_0=1$ so that

$$a_1(n) = \Gamma(1+n) = n!\tag{e1}$$

as expected.

$k=2$

We have $a_2(n) = 2^{n/2} \Gamma \left(\frac{n}{2}+1\right) (c_1 \cos (\pi n)+c_0+s_1 \sin (\pi n))$.

The initial conditions lead to the equations $1=c_0+c_1, 1 = \sqrt{2} (c_0 - c_1)$ giving $c_0= \frac{1}{2}\left(1+\sqrt{\frac{2}{\pi}}\right),c_1 =\frac{1}{2}\left(1-\sqrt{\frac{2}{\pi}}\right)$.

Finally

Complex $n$:

$$a_2(n) = 2^{\frac{n}{2}-1} \Gamma \left(\frac{n}{2}+1\right) \left(2 s_1 \sin (\pi n)\\ +\sqrt{\frac{2}{\pi }} (1-\cos (\pi n))+(\cos (\pi n)+1)\right)\tag{e2a} $$

Graph:

Integer $n$:

For integer we have the simplification $n$ $e^{i \pi n}=(-1)^n$.

Since for even $n$: $\frac{1}{2}\left(\left(1+\sqrt{\frac{2}{\pi }}\right)+\left(1-\sqrt{\frac{2}{\pi }}\right) (-1)^n\right) = 1$ for even $n$ and $=\sqrt{\frac{2}{\pi }}$ for odd $n$ we can also write

$$a_2(n) = 2^{\frac{n}{2}}\Gamma \left(\frac{n}{2}+1\right)(\sqrt{\frac{2}{\pi }})^{\frac{1-(-1)^n}{2}}\tag{e2b} $$

A similar form was given in the OP.

The first few values are

$$\{a_2(n) \}_{n=0}^{n=11} = \{1,1,2,3,8,15,48,105,384,945,3840\}$$

This sequence is listed here https://oeis.org/A006882 as double factorials n!!: a(n) = n*a(n-2) for n > 1, a(0) = a(1) = 1.

$k=3$

This case is even more interesting since $n$ has three values $\mod(3)$

We have, to begin with,

$$a_{3}(n) = 3^{n/3} \left(c_0+c_1 e^{\frac{2 i \pi n}{3}}+c_2 e^{\frac{4 i \pi n}{3}}\right) \Gamma \left(\frac{n}{3}+1\right)\tag{e3a}$$

The inital conditions lead to these linear equations for the coefficients

$$ \begin{array}{c} 1=c_0+c_1+c_2 \\ 1=\sqrt[3]{3} \left(c_0+e^{\frac{2 i \pi }{3}} c_1+e^{-\frac{1}{3} (2 i \pi )} c_2\right) \Gamma \left(\frac{4}{3}\right) \\ 1=3^{2/3} \left(c_0+e^{-\frac{1}{3} (2 i \pi )} c_1+e^{\frac{2 i \pi }{3}} c_2\right) \Gamma \left(\frac{5}{3}\right) \\ \end{array} $$

Real $n$:

The solution of which gives rise to

$$a_3(n) = 3^{n/3} \Gamma \left(\frac{n}{3}+1\right)\left( \frac{1}{3}\left(1+2 (-1)^n \cos \left(\frac{\pi n}{3}\right)\right) \\ +\frac{1}{3} \left(1+(-1)^{2/3} e^{\frac{2 i \pi n}{3}}-\sqrt[3]{-1} e^{\frac{4 i \pi n}{3}}\right) \left(\frac{3^{5/6} \Gamma \left(\frac{1}{3}\right)}{ 4 \pi }\right) \\ +\frac{1}{3}\left((-1)^{n+1} \left(\sqrt{3} \sin \left(\frac{\pi n}{3}\right)+\cos \left(\frac{\pi n}{3}\right)\right)+1\right)\left(\frac{6 \sqrt[6]{3} \Gamma \left(\frac{2}{3}\right)}{ 4 \pi }\right)\right)\tag{e3b} $$

Graph:

Integer $n$:

Simplifying according to the remainder$\mod(3)$ we have

$a_{3}(n) = \left\{ \begin{array}{ll} 3^j \Gamma (j+1) & n=3j + 0 \; \;\;\{\}_{j=0}^{j=5}= \{1,3,18,162,1944,29160\},\; \text{https://oeis.org/A032031}\\ \frac{3^{j+\frac{3}{2}} \Gamma \left(\frac{2}{3}\right) \Gamma \left(j+\frac{4}{3}\right)}{2 \pi } & n=3j+1 \;\;\; \{\}_{j=0}^{j=5} = \{1,4,28,280,3640,58240\},\; \text{https://oeis.org/A007559}\\ \frac{3^{j+\frac{3}{2}} \Gamma \left(\frac{1}{3}\right) \Gamma \left(j+\frac{5}{3}\right)}{4 \pi } & n=3j+2 \;\; \; \{\}_{j=0}^{j=5}=\{1,5,40,440,6160,104720\}, \text{https://oeis.org/A034000}\\ \end{array} \right.\tag{e3c} $

$k=4$

Real $n$:

$$a_{4}(n) = \frac{2^{\frac{n}{2}-2}}{3\pi} \Gamma \left(\frac{n}{4}+1\right) \left(-3 \sqrt{\pi } \left(-1+e^{\frac{i \pi n}{2}}\right) \left(1+e^{i \pi n}\right)+3 \pi \left(1+e^{\frac{i \pi n}{2}}\right) \left(1+e^{i \pi n}\right)-\left(-1+e^{i \pi n}\right) \left(4 i e^{\frac{i \pi n}{2}} \left(\Gamma \left(\frac{5}{4}\right)-2 \Gamma \left(\frac{7}{4}\right)\right)+\Gamma \left(\frac{1}{4}\right)+8 \Gamma \left(\frac{7}{4}\right)\right)\right)$$

Graph:

Integer $n$:

$a_{4}(n) = \left\{ \begin{array}{ll} 4^j \Gamma (j+1) & n=4j + 0 \; \;\;\{\}_{j=0}^{j=5}= \{1,4,32,384,6144,122880\},\; \text{https://oeis.org/A047053}\\ \frac{2^{2 j+\frac{7}{2}} \Gamma \left(\frac{7}{4}\right) \Gamma \left(j+\frac{5}{4}\right)}{3 \pi } & n=4j+1 \;\;\; \{\}_{j=0}^{j=5} = \{1,5,45,585,9945,208845\},\; \text{https://oeis.org/A007696}\\ \frac{2^{2 j+1} \Gamma \left(j+\frac{3}{2}\right)}{\sqrt{\pi }} & n=4j+2 \;\; \; \{\}_{j=0}^{j=5}=\{1,6,60,840,15120,332640\}, \text{https://oeis.org/A000407}\\ \frac{2^{2 j+\frac{3}{2}} \Gamma \left(\frac{1}{4}\right) \Gamma \left(j+\frac{7}{4}\right)}{3 \pi } & n=4j+3 \;\; \; \{\}_{j=0}^{j=5}=\{1, 7, 77, 1155, 21945, 504735\}, \text{https://oeis.org/A034176}\\ \end{array} \right. $

$k=5$

Real $n$: lengthy expression, dropped here

Graph:

Integer $n$:

$a_{5}(n) = \left\{ \begin{array}{ll} 5^j \Gamma (j+1) & n=5j + 0 \; \;\;\{\}_{j=0}^{j=5}= \{1,5,50,750,15000,375000\},\; \text{https://oeis.org/A052562} \\ \frac{\sqrt[10]{\frac{1}{2} \left(25-11 \sqrt{5}\right)} \;5^{j+\frac{11}{5}} \Gamma \left(\frac{9}{5}\right) \Gamma \left(j+\frac{6}{5}\right)}{8 \pi } & n=5j+1 \;\;\; \{\}_{j=0}^{j=5} = \{1,6,66,1056,22176,576576\},\; \text{https://oeis.org/A008548} \\ \frac{\sqrt[10]{\frac{1}{2} \left(25+11 \sqrt{5}\right)}\; 5^{j+\frac{11}{5}} \Gamma \left(\frac{8}{5}\right) \Gamma \left(j+\frac{7}{5}\right)}{12 \pi } & n=5j+2 \;\; \; \{\}_{j=0}^{j=5}=\{1,7,84,1428,31416,848232\}, \text{https://oeis.org/A034323} \\ \frac{5^j \Gamma \left(j+\frac{8}{5}\right)}{\Gamma \left(\frac{8}{5}\right)} & n=5j+3 \;\; \; \{\}_{j=0}^{j=5}=\{1,8,104,1872,43056,1205568\}, \text{https://oeis.org/A034300} \\ \frac{\sqrt[10]{\frac{1}{2} \left(25-11 \sqrt{5}\right)}\; 5^{j+\frac{16}{5}} \Gamma \left(\frac{11}{5}\right) \Gamma \left(j+\frac{9}{5}\right)}{48 \pi } & n=5j+4 \;\; \; \{\}_{j=0}^{j=5}=\{1,9,126,2394,57456,1666224\}, \text{https://oeis.org/A034301} \\ \end{array} \right. $

Discussion

1. Simpler expressions for some $a_{k}(n)$ in terms of the Pochhammer symbol have been given in OEIS by G. C. Greubel, Aug 23 2019.

2. The expressions given here lead somewhat miraculously to integer numbers. The miracle can be solved by using the basic formula for the Gamma function $z \Gamma(z) = \Gamma(z+1)$ and the relation $\Gamma (z) \Gamma (1-z)=\pi \csc (\pi z)$.

3. I compared the results of the two solutions given up to now. For $k=3$, I found that they differ. That of pregunton (https://math.stackexchange.com/a/3488935/198592) agrees with OEIS A007661 "Triple factorial numbers a(n) = n!!!, defined by a(n) = n*a(n-3), a(0) = a(1) = 1, a(2) = 2. Sometimes written n!3." I call this the standard solution. Mine is different and not listed in OEIS. The reason is the different initial conditions. I had assumed that $a_3(2) = 1$, the "standard" is $a_3(2)=2$. For $k=2$ the standard is $a_2(0)=a_2(1)=1$, and consequently my series coincides with the standard. I don't quite understand at the moment why the generalization of the cases $k=1$ and $k=2$ should be a deviation from the constant initial conditions for all values below $k$ (my conditions $(2)$). And what which conditions would be "natural" for $k=4$ and higher $k$? I have found the answer here http://mathworld.wolfram.com/Multifactorial.html. The rule for the initial conditions is simply $a_{m}= m$ for $m=1..k+1$.