I would like to prove the following statement which I got from the same method Euler used to "solve" the Basel problem (using the infinite product of the sine function): $$\sum_{a_1=1}^\infty \ \sum_{a_2=a_1+1}^\infty ... \sum_{a_n=a_{n-1}+1}^\infty \Big(\frac{1}{(a_1a_2... a_n)^2}\Big) = \frac{\pi^{2n}}{(2n+1)!}$$
For $n=1$ it is quite clear, as we end up with $$\sum_{a_1=1}^\infty \Big(\frac{1}{a_1^2}\Big) = \frac{\pi^{2}}{6}$$ which is a result we see pretty often. More interestingly, for $n=2$ we have
$$\sum_{a_1=1}^\infty \ \sum_{a_2=a_1+1}^\infty \Big(\frac{1}{(a_1a_2)^2}\Big) = \frac{\pi^{4}}{120}$$
This nested sum is quite fun and seems to be related to the Zeta function, for example we can quickly find that $\zeta(4)$ is a very simple combination of this sum: $$\zeta(4) = \Big(\sum_{a_1=1}^\infty \Big(\frac{1}{a_1^2}\Big)\Big)^2-2\Big(\sum_{a_1=1}^\infty \ \sum_{a_2=a_1+1}^\infty \Big(\frac{1}{(a_1a_2)^2}\Big)\Big)= \Big(\frac{\pi^{2}}{6}\Big)^2-2\times \frac{\pi^{4}}{120} = \frac{\pi^4}{90}$$
However, trying to attempt anything with $n>2$ makes everything extremely difficult, and I have absolutely no intuition for how to prove the initial statement
I am still a student in a French secondary school so I'm certain I don't know most of the tricks used in sums. Sorry for any bad grammar and thank you very much for any help!
The key idea is that Newton's identities also work with an infinite number of variables, if all the sums converge. The proof of Newton's identities can easily be generalized to an infinite number of variables if everything is "well-behaved". Let $$ x_j = \frac{1}{j^2} $$ You want to show $$ e_n(x_1,x_2,x_3,\ldots) = \sum _{j_1<j_2<\ldots < j_n} x_{j_1} x_{j_2}\cdots x_{j_n} = \frac{\pi^{2n}}{(2n+1)!} $$ where $e_n(x_1,x_2,x_3,\ldots)$ is the $n$-th elementary symmetric polynomial of $x_1,x_2,\ldots$ (generalized for an infinite number of variables).
Likewise, we have $$ p_k(x_1,x_2,\ldots) = \sum_{j=1}^{\infty}x_j^k = \sum_{j=1}^{\infty} \frac{1}{j^{2k}} = \zeta(2k) $$
Newton's identities state that $$ ne_n(x_1,x_2,\ldots) = \sum_{k=1}^{n}(-1)^{k-1}e_{n-k}(x_1,x_2,\ldots)p_k(x_1,x_2,\ldots) $$ We must check if this equation holds if we plug in the known values for $p_i$ and the assumed values for $e_i.$ We can use this equation for a proof by induction, then. We already know that your claim holds for $n=1$, and we can make the induction step from $n-1$ to $n$ with Newton's identities, if they hold for our assignments of the $e_i$ and the $p_i.$
Well then, we have to show the following $$ \frac{n\pi^{2n}}{(2n+1)!} = \sum_{k=1}^{n}(-1)^{k-1} \frac{\pi^{2n-2k}}{(2n-2k+1)!}\,\zeta(2k) $$ I have not found a simple and elegant way to show this. The following is quite tedious and uses many formulae containing Bernoulli numbers. However, in the end it turns out that your claim is correct.
Divide both sides by $\pi^{2n}$ to get $$ \frac{n}{(2n+1)!} = \sum_{k=1}^{n} \frac{(-1)^{k-1}}{(2n-2k+1)!} \cdot\frac{\zeta(2k)}{\pi^{2k}} $$ Replace the zeta function with Bernoulli numbers: $$ \frac{n}{(2n+1)!} = \sum_{k=1}^{n} \frac{(-1)^{k-1}}{(2n-2k+1)!} \cdot\frac{(-1)^{k+1}(2\pi)^{2k}B_{2k}}{2((2k)!)\pi^{2k}} $$ Multiply with $(2n+1)!$ and simplify $$ n = \sum_{k=1}^{n} \binom{2n+1}{2k} 2^{2k-1}B_{2k} $$ Use the fact that $B_3=B_5 = \ldots = 0$ and perform re-indexing. Multiply with $2$ : $$ 2n = \sum_{k=2}^{2n+1} \binom{2n+1}{k} 2^k B_k $$ Replace $2^{k}$ with $\left(\frac{1}{2}\right)^{-k}$ and multiply with $\left(\frac{1}{2}\right)^{2n+1}$ : $$ \left(\frac{1}{2}\right)^{2n} n = \sum_{k=2}^{2n+1} \binom{2n+1}{k} \left(\frac{1}{2}\right)^{2n+1-k} B_k $$ Add the terms for $k=0$ and $k=1$ $$ \left(\frac{1}{2}\right)^{2n+1} B_0 +\left(\frac{1}{2}\right)^{2n} (2n+1) B_1 + \left(\frac{1}{2}\right)^{2n} n = \sum_{k=0}^{2n+1} \binom{2n+1}{k} \left(\frac{1}{2}\right)^{2n+1-k} B_k $$ With the Bernoulli polynomial $B_{2n+1}(x),$ we get $$ \left(\frac{1}{2}\right)^{2n+1} B_0 +\left(\frac{1}{2}\right)^{2n} (2n+1) B_1 + \left(\frac{1}{2}\right)^{2n} n = B_{2n+1}\left(\frac{1}{2}\right) $$ The right side is $0$, because $B_m\left(\frac{1}{2}\right)$ is a multiple of $B_m,$ which is zero for odd $m >2.$ So the only thing left to show is $$ \left(\frac{1}{2}\right)^{2n+1} B_0 +\left(\frac{1}{2}\right)^{2n} (2n+1) B_1 + \left(\frac{1}{2}\right)^{2n} n = 0 $$ Multiply with $2^{2n+1}$ : $$ B_0 +2 (2n+1) B_1 + 2 n = 0 $$ which is trivial, because $B_0 =1$ and $B_1 = - \frac{1}{2}.$