In example 4.18 of the book Introduction to Probability Models by Sheldon Ross (10th edition), he shows that the probability a symmetric 1-d random talk will come back to the origin in $2n$ steps to be:
$$P^{2n}_{00} = {2n \choose n}\left(\frac{1}{4}\right)^n$$
And for the 2-d case, he goes on to prove:
$$P^{2n}_{00} = {2n \choose n}{2n \choose n}\left(\frac{1}{4}\right)^{2n}$$
This suggests the tantalizing possibility that for the 3-d case we must have:
$$P^{2n}_{00} = {2n \choose n}{2n \choose n}{2n \choose n}\left(\frac{1}{6}\right)^{2n}$$
And in general for $k$ dimensions:
$$P^{2n}_{00} = {2n \choose n}^k\left(\frac{1}{2k}\right)^{2n}$$
I tried to prove this for the 3-d case using the same strategy he employed for the results involving the 1-d and 2-d cases. But with two floating indices, the algebra became quite difficult. This is probably why he omitted it in the book. I'm wondering if there is some way to prove (or disprove) the third equation.
My attempt
Similar to what he did for the 1-d and 2-d cases, we imagine $i$ steps along the positive x-axis and the same amount along the negative $x$. Similarly $j$ along the positive and negative y-axis and $k$ along the positive and negative z-axis. We get:
$$P^{2n}_{00} =\sum_{i,j,k: i+j+k=n} \frac{(2i+2j+2k)!}{(i!)^2 (j!)^2 (k!)^2}\left(\frac{1}{6}\right)^{2n}$$
Substituting $i+j+k=n$ we get:
$$P^{2n}_{00} =\sum_{i,j: i+j\leq n} \frac{(2n)!}{(i!)^2 (j!)^2 ((n-i-j)!)^2}\left(\frac{1}{6}\right)^{2n}$$
It isn't clear how to handle this double summation from here.
The 2D argument does not generalize.
Instead of taking steps in 2D that are $(+1,+0)$, $(-1,+0)$, $(+0,+1)$, and $(+0,-1)$, we can rotate the grid 45 degrees and make steps of $(+1,+1)$, $(+1,-1)$, $(-1,+1)$, and $(-1,-1)$. With this setup, it's easy to see that the changes in $x$-coordinate are independent random variables from the the changes in $y$-coordinate, so we can simply square the 1D probability.
Similarly, $\binom{2n}n^3 8^{-2n}$ is the return probability for the 3D random walk where every step has $8$ possibilities: to change each coordinate simultaneously by $+1$ or $-1$. However, this is no longer equivalent to the random walk where we randomly select one of the coordinates to change, so it doesn't help us! The coincidence that helped us in 2D does not recur in any higher dimension.
Actually, there is a further problem with the proposed formula: with an exponential term of $6^{-2n}$, the formula gives values larger than $1$ for large values of $n$, because it's counterbalanced by $\binom{2n}{n}^3$ or roughly $2^{6n}$.