Introduction
In Wolfram Functions there is a relation between the Riemann Zeta Function $\zeta\left( z \right)$ and the Generalized Hypergeometric Function $\operatorname{_{p}F_{q}}\left( a_{1},\, \dots,\, a_{p};\, b_{1},\, \dots,\, b_{q};\, z \right)$ given by this formulas:
$$
\begin{align*}
\zeta\left( n \right) &= \frac{1}{1 – 2^{1 – n}} \cdot \operatorname{_{n + 1}F_{n}}\left( 1,\, a_{1},\, a_{2},\, \dots,\, a_{n};\, a_{1} + 1,\, a_{2} + 1,\, \dots,\, a_{n} + 1;\, -1 \right)\, /;&\, a_{1} = a_{2} = \dots = a_{n} = 1 \wedge n – 1 \in \mathbb{N}\\
\zeta\left( n \right) &= \frac{1}{1 – 2^{1 – n}} \cdot \operatorname{_{n + 1}F_{n}}\left( \begin{matrix} 1 & a_{1} & a_{2} & \dots & a_{n}\\ a_{1} + 1 & a_{2} + 1 & a_{3} + 1 & \dots & a_{n} + 1 \end{matrix};\, -1 \right)\, /;&\, a_{1} = a_{2} = \dots = a_{n} = 1 \wedge n – 1 \in \mathbb{N} \tag{1.}\\
\end{align*}
$$
and
$$
\begin{align*}
\zeta\left( n \right) &= \operatorname{_{n + 1}F_{n}}\left( 1,\, a_{1},\, a_{2},\, \dots,\, a_{n};\, a_{1} + 1,\, a_{2} + 1,\, \dots,\, a_{n} + 1;\, 1 \right)\, /;&\, a_{1} = a_{2} = \dots = a_{n} = 1 \wedge n – 1 \in \mathbb{N}\\
\zeta\left( n \right) &= \operatorname{_{n + 1}F_{n}}\left( \begin{matrix} 1 & a_{1} & a_{2} & \dots & a_{n}\\ a_{1} + 1 & a_{2} + 1 & a_{3} + 1 & \dots & a_{n} + 1 \end{matrix};\, 1 \right)\, /;&\, a_{1} = a_{2} = \dots = a_{n} = 1 \wedge n – 1 \in \mathbb{N} \tag{2.}\\
\end{align*}
$$
Now I'm wondering how they came up with that, since it doesn't exactly seem trivial to me. I've never seen these relations before and Wolfram Functions doesn't give any references, hence the question: How can the relation of the Riemann Zeta Function to the Generalized Hypergeometric Function be proven / derived?
My First Idea
The formula $1.$ reminds me of one of a series expansions for the Riemann Zeta Function given by this formula:
$$
\begin{align*}
\zeta\left( x \right) &= \frac{1}{1 – 2^{1 – x}} \cdot \sum_{k = 1}^{\infty} \left( \frac{\left( -1 \right)^{k – 1}}{k^{x}} \right),\, \text{ for } 0 < \Re\left( x \right) \ne 1\\
\end{align*}
$$
With some manipulation we get $\frac{1}{1 – 2^{1 – n}} \cdot \operatorname{_{n + 1}F_{n}}\left( \begin{matrix} 1 & a_{1} & a_{2} & \dots & a_{n}\\ a_{1} + 1 & a_{2} + 1 & a_{3} + 1 & \dots & a_{n} + 1 \end{matrix};\, -1 \right) = \frac{1}{1 – 2^{1 – n}} \cdot \sum_{k = 1}^{\infty} \left( \frac{\left( -1 \right)^{k – 1}}{k^{n}} \right)$ and $\operatorname{_{n + 1}F_{n}}\left( \begin{matrix} 1 & a_{1} & a_{2} & \dots & a_{n}\\ a_{1} + 1 & a_{2} + 1 & a_{3} + 1 & \dots & a_{n} + 1 \end{matrix};\, -1 \right) = \sum_{k = 1}^{\infty} \left( \frac{\left( -1 \right)^{k – 1}}{k^{n}} \right)$.
Via using the formula $\operatorname{_{p}F_{q}}\left( a_{1},\, \dots,\, a_{p};\, b_{1},\, \dots,\, b_{q};\, z \right) = \sum_{k = 0}^{\infty} \left( \frac{\left( a_{1} \right)_{k} \cdot \left( a_{2} \right)_{k} \cdots \left( a_{p} \right)_{k}}{\left( b_{1} \right)_{k} \cdot \left( b_{2} \right)_{k} \cdots \left( b_{q} \right)_{k}} \cdot \frac{z^{k}}{k!} \right)$ where $\left( x \right)_{y}$ Pochhammer Symbol or Rising Factorial ($\left( x \right)_{y} \equiv \Gamma\left( x + y \right) / \Gamma\left( x \right)$) we get:
$$
\begin{align*}
\sum_{k = 0}^{\infty} \left( \frac{\left( a_{1} \right)_{k} \cdot \left( a_{2} \right)_{k} \cdots \left( a_{n + 1} \right)_{k}}{\left( a_{1} + 1 \right)_{k} \cdot \left( a_{2} + 1 \right)_{k} \cdots \left( a_{n} + 1 \right)_{k}} \cdot \frac{\left(
-1 \right)^{k}}{k!} \right) &= \sum_{k = 1}^{\infty} \left( \frac{\left( -1 \right)^{k – 1}}{k^{n}} \right)\\
\sum_{k = 1}^{\infty} \left( \frac{\left( a_{1} \right)_{k} \cdot \left( a_{2} \right)_{k} \cdots \left( a_{n + 1} \right)_{k}}{\left( a_{1} + 1 \right)_{k} \cdot \left( a_{2} + 1 \right)_{k} \cdots \left( a_{n} + 1 \right)_{k}} \cdot \frac{\left(
-1 \right)^{k}}{k!} \right) &= \sum_{k = 1}^{\infty} \left( \frac{\left( -1 \right)^{k – 1}}{k^{n}} \right) + \frac{\left( a_{1} \right)_{1} \cdot \left( a_{2} \right)_{1} \cdots \left( a_{n + 1} \right)_{1}}{\left( a_{1} + 1 \right)_{1} \cdot \left( a_{2} + 1 \right)_{1} \cdots \left( a_{n} + 1 \right)_{1}}\\
\end{align*}
$$
But I don't see any real transformation from here that could help here.
My Second Idea
My second idea would be to look at the equations somehow using differential equations. Equations $\operatorname{_{n + 1}F_{n}}\left( a_{1},\, \dots,\, a_{n + 1};\, b_{1},\, \dots,\, b_{n};\, z \right)$ and $z^{1 – b_{1}} \cdot \operatorname{_{n + 1}F_{n}}\left( 1 + a_{1} – b_{1},\, 1 – a_{2} – b_{1},\, \dots,\, 1 + a_{n + 1} – b_{1};\, 2 – b_{1},\, 1 – b_{2} – b_{1},\,\dots,\, 1 – b_{n} – b_{1};\, z \right)$ are two known solutions of the differential equation $\left[ \vartheta\left( \vartheta + b_{1} – 1 \right) \cdots \left( \vartheta + b_{n} – 1 \right) – z\left( \vartheta + a_{1} \right) \cdots \left( \vartheta + a_{n + 1} \right) \right] y = 0$ (Bailey 1935, p. 8), where $\vartheta = z \cdot \frac{\operatorname{d}}{\operatorname{d}z}$ (Rainville 1971, Koepf 1998, p. 27).
If we were to linearize the ODEs $\left[ \vartheta\left( \vartheta + b_{1} – 1 \right) \cdots \left( \vartheta + b_{n} – 1 \right) \pm \left( \vartheta + a_{1} \right) \cdots \left( \vartheta + a_{n + 1} \right) \right] y = 0$, we could use the Powerseries Method to solve it, with which we might find a solution using the Riemann zeta function but on the one hand I think that linearizing would be a bit exhausting, but also that I don't know enough special series developments to manipulate it that way.
My Third Idea
The series representation of $\zeta\left( z \right)$ at $z = 1$ is given by $\frac{1}{z – 1} + \sum_{k = 1}^{\infty} \left( \frac{\left( -1 \right)^{k}}{k!} \cdot \gamma_{k} \cdot \left( z – 1 \right)^{k} \right)$ where $\gamma_{k}$ is the Stieltjes Constants aka
$$
\begin{align*}
\zeta\left( n + 1 \right) &= \frac{1}{n} + \sum_{k = 1}^{\infty} \left( \frac{\left( -1 \right)^{k}}{k!} \cdot \gamma_{k} \cdot n^{k} \right).
\end{align*}
$$
If we now substitute $c_{k} := \left( -1 \right)^{k} \cdot \gamma_{k}$ aka $\zeta\left( n + 1 \right) = \frac{1}{z} + \sum_{k = 1}^{\infty} \left( c_{k} \cdot \frac{n^{k}}{k!} \right)$, then we can apply the formula $\frac{c_{k + 1}}{c_{k}} = \frac{\left( k + a_{1} \right) \cdot \left( k + a_{2} \right) \cdots \left( k + a_{p} \right)}{\left( k + b_{1} \right) \cdot \left( k + b_{2} \right) \cdots \left( k + b_{q} \right) \cdot \left( k + 1 \right)}$ to get:
$$
\begin{align*}
\frac{c_{k + 1}}{c_{k}} &= \frac{\left( -1 \right)^{k + 1} \cdot \gamma_{k + 1}}{\left( -1 \right)^{k} \cdot \gamma_{k}}\\
\frac{c_{k + 1}}{c_{k}} &= -\frac{\left( -1 \right)^{k} \cdot \gamma_{k + 1}}{\left( -1 \right)^{k} \cdot \gamma_{k}}\\
\frac{c_{k + 1}}{c_{k}} &= -\frac{\gamma_{k + 1}}{\gamma_{k}}\\
\end{align*}
$$
I don't see any further logical manipulation here.
Best Answer
For the first formula, we need to prove
$$\zeta(n)=\sum\limits_{k=0}^{\infty} \frac{\left(\prod_{j=1}^{n+1} (1)_k\right)}{\left(\prod _{j=1}^n (2)_k\right)} \frac{1^k}{k!}\tag{1}$$
where the right-side is the expansion of the first hypergeometric $\operatorname{_{p}F_{q}}$ function.
Noting that
$$\frac{(1)_k}{(2)_k}=\frac{1}{k+1}\tag{2}$$
and
$$(1)_k=k!\tag{3}$$
leads to
$$\frac{\prod _{j=1}^{n+1} (1)_k}{\prod _{j=1}^n (2)_k}=\frac{k!}{(k+1)^n}\tag{4}$$
so formula (1) above simplifies to
$$\zeta(n)=\sum\limits_{k=0}^{\infty} \frac{k!}{(k+1)^n} \frac{1^k}{k!}=\sum\limits_{k=1}^{\infty} \frac{1}{k^n}\tag{5}$$
For the second formula, we need to prove
$$\zeta(n)=\frac{1}{1-2^{1-n}} \sum\limits_{k=0}^{\infty} \frac{\left(\prod_{j=1}^{n+1} (1)_k\right)}{\left(\prod _{j=1}^n (2)_k\right)} \frac{(-1)^k}{k!}\tag{6}$$
where the right-side is the expansion of the second hypergeometric $\operatorname{_{p}F_{q}}$ function.
Using the relationship in formula (4) above formula (6) above simplifies to
$$\zeta(n)=\frac{1}{1-2^{1-n}} \sum\limits_{k=0}^{\infty} \frac{k!}{(k+1)^n} \frac{(-1)^k}{k!}=\frac{1}{1-2^{1-n}} \sum\limits_{k=1}^{\infty} \frac{(-1)^{k+1}}{k^n}\tag{7}$$