Integral involving hypergeometric function $\int_0^1[{}_2F_1(\frac13,\frac23;1;x^3)]^2dx$

Question: How to prove $$I=\int_0^1\bigg[{}_2F_1\left(\frac13,\frac23;1;x^3\right)\bigg]^2dx=\frac{\sqrt3}{32\pi^5}\Gamma\left(\frac13\right)^9?$$

Source: An integral competition post of my country.

Attempt
Recall the series definition of hypergeometric function $$_2F_1(a,b,c,x)=\sum_{n=0}^\infty\frac{(a)_n(b)_n}{(c)_nn!}x^n,$$ we can transform $I$ into the series form $$I=\sum_{n,m=0}^\infty\frac{(a)_n(a)_m(b)_n(b)_m}{(c)_n(c)_mn!m!(3n+3m+1)}$$ But I can not handle this series.
I also thought of using complex method. $$I=\int_0^1\bigg[{}_2F_1\left(\frac13,\frac23;1;x\right)\bigg]^2\frac{x^{-2/3}}3dx$$ then let $f(z)=\bigg[{}_2F_1\left(\frac13,\frac23;1;z\right)\bigg]^2(-z)^{-2/3}$ and use keyhole contour, where $(\cdot)^{-2/3}$ is the principal branch of the multi-valued function. But the nature of the branch of the integrand in the inteval $[1,\infty)$ is too complex for me to handle. It involves another definite integral which is similar to $I$.

Before attacking the integral, I mention something about cubic theta function. The whole solution heavily exploits tools from modular forms. The "footnote" contains more information.

The three cubic theta functions are defined by $$\begin{aligned} a(q) &= \sum_{m,n} q^{m^2+mn+n^2}\\ b(q) &= \sum_{m,n} \zeta_3^{m-n} q^{m^2+mn+n^2}\\ c(q) &= \sum_{m,n} q^{{(m+\frac{1}{3})^2+(m+\frac{1}{3})(n+\frac{1}{3})+(n+\frac{1}{3})^2}} \end{aligned}$$ where $\zeta_3 = e^{2\pi i/3}$, sum is over all $m,n\in \mathbb{Z}$. Then it can be shown$^1$ that $$a(q)^3 = b(q)^3+c(q)^3$$ $$a(q) = \frac{\eta^3(q) + 9 \eta^3(q^9)}{\eta (q^3)}\qquad b(q) = \frac{\eta^3(q)}{\eta(q^3)}\qquad c(q) = 3\frac{\eta^3(q^3)}{\eta(q)}$$ where $\eta(q) = q^{1/24} \prod_{n\geq 1}(1-q^n)$ is the Dedekind eta function.

Define $$K_3(m) = {_2F_1}(\frac{1}{3},\frac{2}{3};1;m) $$ Similar to elliptic integrals, denote $K_3'(m) = K_3(1-m), m' = 1-m$. Then one easily shows (I omit the subscript $3$): $$\frac{d}{dm}(\frac{K'}{K}) = -\frac{\sqrt{3}}{2\pi}\frac{1}{mm'K^2}$$

Moreover, letting $q= \exp(-\frac{2\pi}{\sqrt{3}}\frac{K'(m)}{K(m)})$, the following inversion formula holds$^2$ when $0<m<1$: $$a(q) = K(m)\qquad b(q)=(1-m)^{1/3} K(m)\qquad c(q) = m^{1/3} K(m)$$

Now we tackle the integral, $$I = \frac{1}{3}\int_0^1 {{m^{ - 2/3}}K{{(m)}^2}dm} $$ we make the substitution $q = \exp ( - \frac{{2\pi }}{{\sqrt 3 }}\frac{{K'(m)}}{{K(m)}})$, the above formulas imply $dq = \frac{{q}}{{mm'{K^2}}}dm$, as $m$ increases from $0$ to $1$, $q$ increases from $0$ to $1$. $$I = \frac{1}{3}\int_0^1 {\frac{{b{{(q)}^3}c(q)}}{{mm'{K^2}}}dm} = \frac{1}{3}\int_0^1 {\frac{{b{{(q)}^3}c(q)}}{q}dq} = \int_0^1 {\frac{{\eta {{(q)}^8}}}{q}dq} $$ Next, I will use notation $\eta(q),\eta(\tau)$ interchangably (the common notation in context of modular forms), where $q = e^{2\pi i \tau}$. Make $q=e^{-2\pi x}$, then $I$ becomes $$I = 2\pi \int_0^\infty {\eta {{(ix)}^8}dx} = 2\pi \int_0^\infty {{x^2}\eta {{(ix)}^8}dx} $$ where in last step, I used $\eta(-1/\tau) = \sqrt{-i\tau} \eta(\tau)$. Transform it back to $q$: $$\tag{1} I = \frac{1}{{4{\pi ^2}}}\int_0^1 {\frac{{{{\ln }^2}q}}{q}\eta {{(q)}^8}dq} $$ It can be shown that$^3$: $$\eta {(q)^8} = - \frac{1}{2}\sum\limits_{v \in S} {({v_0} - {v_1})({v_1} - {v_2})({v_0} - {v_2}){q^{{{\left\| v \right\|}^2}/6}}}$$ $$S = \left\{ {v \in {\mathbb{R}^3}|v = ({v_0},{v_1},{v_2}) = (3n,3m + 1,3r - 1),n + m + r = 0,n,m,r\in\mathbb{Z}} \right\}$$ with $\|v\|$ the norm of a vector. Plug this into (1): $$I = \frac{{ - 1}}{{{{(2\pi )}^2}}}{6^3}\sum\limits_{v \in S} {\frac{{({v_0} - {v_1})({v_1} - {v_2})({v_0} - {v_2})}}{{{{\left\| v \right\|}^6}}}} $$ Denote $\rho = e^{\pi i/3}$. Note that $({v_0} - {v_1})({v_1} - {v_2})({v_0} - {v_2}) = 2\Re {({v_0} + \rho {v_1})^3}$ and $${\left\| v \right\|^6} = 8{({v_0} + \rho {v_1})^3}{({v_0} + {\rho ^{ - 1}}{v_1})^3}$$ we obtain $$I = \frac{{ - 27}}{{2{\pi ^2}}}\Re \sum\limits_{v\in S} {\frac{1}{{{{({v_0} + {\rho ^{ - 1}}{v_1})}^3}}}} = - \frac{{27}}{{2{\pi ^2}}}\Re \sum\limits_{(m,n) \in {\mathbb{Z}^2}} {\frac{1}{{{{(3n + {\rho ^{ - 1}}(3m + 1))}^3}}}}$$ The latter can be recognized as an Eisenstein series of level $3$, but to calculate its value, it is best to use Weierstrass elliptic function. Let $\wp_{1,\rho}$ denote this elliptic function with periods $\{1,\rho\}$, then $${\wp _{1,\rho }}'(z) = - 2\sum\limits_{n,m} {\frac{1}{{{{(z + n + m\rho )}^3}}}} $$ gives $$I=\frac{1}{{4{\pi ^2}}}\Re \left[{\wp _{1,\rho }}'(\frac{{{\rho ^{ - 1}}}}{3})\right] = \frac{{{\omega ^3}}}{{4{\pi ^2}}}\Re\left[ {\wp _{\omega ,\omega \rho }}'(\frac{{{\omega\rho ^{ - 1}}}}{3})\right]$$ where $\omega = \Gamma(1/3)^3/(2\pi)$, then it is well-known that modular invariants associated to periods $\{\omega,\omega\rho\}$ are $g_2 = 0, g_3 = 1$. Therefore ${\wp _{\omega ,\omega \rho }}'(\frac{{\omega {\rho ^{ - 1}}}}{3})$ is the $y$-coordinate of a $3$-torsion of the elliptic curve $y^2 = 4x^3 - g_2 x - g_3 = 4x^3 -1$, which can be readily calculated to be $\sqrt{3}$. Finally we finish the calculation:$I = \omega^3\sqrt 3/(4\pi^2)$.

$^1$: Proof outline: $a(q^3),b(q^3),c(q^3)$ are modular forms of weight $1$ and level $27$, therefore it suffices to verify their $q$-expansions to certain power of $q$. A self-contained approach can be found in the 1994 paper Cubic Analogues of the Jacobian Theta Function.

$^2$: Proof outline: $f=c^3(\tau)/a^3(\tau)$ is modular function of $\Gamma_0(3)$, by a fact in modular forms, $b(\tau)$ satisfies a 2nd order ODE in terms of $f$, its coefficients are rational functions of $f$ since modular curve $X(3)$ has genus $0$. Therefore, in certain region of $\mathbb{H}$, $b(\tau) = (1-f)^{1/3} K_3(f)$, we could replace $\tau$ by $\gamma\tau$ for $\gamma\in \Gamma_0(3)$, modularity of $b$ allows us to isolate the $\tau$. But doing this replacement might change it into another linear independent solution of the ODE, which explains why $K'/K$ arises. The details are more delicate.

$^3$: The exponent $8$ is special here, which is the dimension of semisimple Lie algebra $A_2$. There is a corresponding formula for $\eta(q)^d$ each semisimple Lie algebra with dimension $d$. See Affine Root Systems and Dedekind's eta-Function by I.G. Macdonald.

Not an answer, but an extended comment for now.

This hypergeometric function is a special case, and some complicated quadratic and cubic transformations apply to it. See this like for reference: https://dlmf.nist.gov/15.8.

The formulas 15.8.25 and 15.8.26 both apply here.

However, the most interesting one is so called Ramanujan’s Cubic Transformation (15.8.33):

$${_2 F_1} \left( \frac13, \frac23;1;x^3 \right)= \frac{1}{1+2 x} {_2 F_1} \left( \frac13, \frac23;1;1- \frac{(1-x)^3}{(1+2x)^3}\right)$$

Update:

The iteration:

$$x_{n+1}=\left( 1- \frac{(1-x_n)^3}{(1+2x_n)^3}\right)^{1/3}$$

Converges to $x_{\infty}=1$ for any $x \in (0,1]$. Not sure how to use this, because ${_2 F_1} \left( \frac13, \frac23;1;1 \right)= \infty$.

This transformation is related to the cubic analogue of the arithmetic-geometric mean. See the references at DLMF and also these questions:

Integral identity related with cubic analogue of arithmetic-geometric mean

Evaluate the integral $\int_0^\infty \frac{dx}{\sqrt{(x^3+a^3)(x^3+b^3)}}$

Some formulas from the question above (and Nemo's answer) might be useful here, for example:

$$\int_0^\infty \frac{dt}{\sqrt{(t^3+1)(t^3+p)}}=\frac{2 \pi}{3 \sqrt{3}} {_2F_1} \left(\frac{1}{2},\frac{2}{3};1;1-p \right)= \\ =\frac{2 \pi}{3 \sqrt{3}p^{1/3}}{_2F_1} \left(\frac{1}{3},\frac{2}{3};1;\frac{(1-\sqrt{p})^2}{-4\sqrt{p}} \right)$$

This is just an application of already linked transformations, and can be applied backwards in this case.