There are a couple things that could be slowing this down, but the main one is likely that the Sigma symbol parses into sum() which is much slower than the more straight forward add(). Try defining J this way:
J := unapply(8*Pi^(3/2)*r*R*(add((2*r*R)^(2*i)*pochhammer((1/2)*n, i)
*hypergeom([2*j+2*l+1, .5], [2*j+2*l+1.5], -1)*(.5
*hypergeom([-n-2*i-l+1.5, l+.5], [l+1.5], sin(arctan(-Z/sqrt(R^2+r^2)))^2)/(2*l+1))
l = 0 .. 100))/(factorial(i-j)*factorial(j)),
j = 0 .. i))/factorial(i), i = 0 .. 100)),[n,phi]);
Your worksheet still takes about three minutes to execute, but I killed the other version after 10 minutes.