Prove that $_4F_3\left(\frac13,\frac13,\frac23,\frac23;1,\frac43,\frac43;1\right)=\frac{\Gamma \left(\frac13\right)^6}{36 \pi ^2}$ I found an interesting problem about generalized hypergeometric series in MO, stating that:

$$\, _4F_3\left(\frac{1}{3},\frac{1}{3},\frac{2}{3},\frac{2}{3};1,\frac{4}{3},\frac{4}{3};1\right)=\sum_{n=0}^\infty \left(\frac{(\frac13)_k (\frac23)_k}{(1)_k (\frac43)_k}\right)^2=\frac{\Gamma \left(\frac{1}{3}\right)^6}{36 \pi ^2}$$

This numerically agrees, but I found no proof using either elementary properties of hypergeometric functions (e.g. cyclic sum) or classical Gamma formulas (e.g. Dougall formula). I bet it has something to do with modular forms and elliptic $K$ integral, but the exact relation remain elusive.
How to prove this identity? What will be its motivation? Can we generate other Gamma evaluation of high order hypergeometric series using the method of proving it? Any help will be appreciated.
 A: Wow, amazing! Solved 9 years later! Thanks all for digging this up, and then for solving it. Can this give a general form for
$$_4F_3(\frac1m,\frac1m,\frac2m,\frac2m;\frac{m+1}m,\frac{m+1}m,1;1)$$
I should probably give some motivations for this. In the following paper, I looked at the expected exit time of a planar Brownian motion starting at 0 from a regular $m$-gon centered at 0:
https://projecteuclid.org/euclid.ecp/1465262013
It is (up to a constant which depends on the size of the polygon)
$$_4F_3(\frac1m,\frac1m,\frac2m,\frac2m;\frac{m+1}m,\frac{m+1}m,1;1)\times \frac{m^2}{\beta(1/m,(m-2)/m)^2},$$
which doesn't exactly roll off the tongue. However, for an equilateral triangle there is a different method for calculating this, and it gives $1/6$. So we get an identity by equating the two, and that is the identity. Now, the question is, can we use this method to get a nicer expression for the $_4F_3$ for larger $m$? This would then be a nicer expression for the expected exit time of Brownian motion from the regular $m$-gon.
A purely analytic (i.e. not probabilistic) version of this all can be found here, because the expected exit time is basically the Hardy H^2 norm of the domain, up to a constant.
https://arxiv.org/abs/1205.2458
