Evaluating $\int_0^1\frac{x^{2/3}(1-x)^{-1/3}}{1-x+x^2}dx$

How can we prove $$\int_0^1\frac{x^{2/3}(1-x)^{-1/3}}{1-x+x^2}\mathrm{d} x=\frac{2\pi}{3\sqrt 3}?$$

Thought 1
It cannot be solved by using contour integration directly. If we replace $-1/3$ with $-2/3$ or $1/3$ or something else, we can use contour integration directly to solve it.
Thought 2
I have tried substitution $x=t^3$ and $x=1-t$. None of them worked. But I noticed that the form of $1-x+x^2$ does not change while applying $x=1-t$.
Thought 3
Recall the integral representation of $_2F_1$ function, I was able to convert it into a formula with $_2F_1\left(2/3,1;4/3; e^{\pi i/3}\right)$ involved. But I think it will only make the integral more "complex". Moreover, I prefer a elementary approach. (But I also appreciate hypergeometric approach)

Update

I have now finally found a way to take my hypergeometric solution all the way to its final elementary form.

Let $$I = \int_0^1 \frac{x^{2/3}}{\sqrt[3]{1 - x} (1 - x +x^2)} \, dx. \tag1$$ Enforcing a substitution of $x \mapsto 1 - x$ leads to $$I = \int_0^1 \frac{(1 - x)^{2/3}}{\sqrt[3]{x} (1 - x + x^2)} \, dx. \tag2$$ Adding (1) to (2) produces $$I = \frac{1}{2} \int_0^1 \frac{dx}{\sqrt[3]{x - x^2} (1 - x + x^2)}.$$

Expanding the second term appearing in the denominator in terms of a geometric series, we have \begin{align} I &= \frac{1}{2} \int_0^1 \frac{dx}{(x - x^2)^{1/3} [1 - (x - x^2)]}\\ &= \frac{1}{2} \int_0^1 \frac{1}{(x - x^2)^{1/3}} \sum_{n = 0}^\infty (x - x^2)^n \, dx\\ &= \frac{1}{2} \sum_{n = 0}^\infty \int_0^1 x^{n - 1/3} (1 - x)^{n - 1/3} \, dx\\ &= \frac{1}{2} \sum_{n = 0}^\infty \operatorname{B} \left (n + \frac{2}{3}, n + \frac{2}{3} \right ), \tag3 \end{align} where $\operatorname{B}(x,y)$ is the Beta function. Making use of the result $$\operatorname{B} (x,x) = \frac{\sqrt{\pi} 2^{1 - 2x} \Gamma (x)}{\Gamma \left (x + \frac{1}{2} \right )},$$ the sum in (3) can be written as \begin{align} I &= \frac{\sqrt{\pi}}{2 \sqrt[3]{2}} \sum_{n = 0}^\infty \frac{\Gamma \left (n + \frac{2}{3} \right )}{\Gamma \left (n + \frac{7}{6} \right ) 4^n}\\ &= \frac{\sqrt{\pi}}{2 \sqrt[3]{2}} \cdot \frac{\Gamma (\frac{2}{3})}{\Gamma (\frac{7}{6})} \sum_{n = 0}^\infty \frac{\left (\frac{2}{3} \right )_n (1)_n}{\left (\frac{7}{6} \right )_n 4^n n!}\\ &= \frac{\sqrt{\pi}}{2 \sqrt[3]{2}} \cdot \frac{\Gamma (\frac{2}{3})}{\Gamma (\frac{7}{6})}\ _2F_1 \left (\frac{2}{3}, 1; \frac{7}{6}; \frac{1}{4} \right ),\tag4 \end{align} where $_2F_1 (a,b;c;x)$ is the Gauss hypergeometric function.

To reduce the hypergeometric function that appears in (4) into elementary form, we proceed as follows.

Firstly, since $_2F_1 (a,b;c;x) =\ _2F_1 (b,a;c;x)$ on applying the second of Pfaff's transformations, namely $$_2F_1 (a,b;c;x) = (1 - x)^{-a}\ _2F_1 \left (a,c-b;c;\frac{x}{x - 1} \right ),$$ to the hypergeometric function, we have $$_2F_1 \left (1, \frac{2}{3}; \frac{7}{6}; \frac{1}{4} \right ) = \frac{4}{3}\ _2F_1 \left (1, \frac{1}{2}; \frac{7}{6}; -\frac{1}{3} \right ).\tag5$$ Next, applying Euler's transformation, namely $$_2F_1 (a,b;c;x) = (1 - x)^{c - a - b}\ _2F_1 (c-a,c-b;c;x),$$ we have $$_2F_1 \left (1, \frac{2}{3}; \frac{7}{6}; \frac{1}{4} \right ) = \frac{4^{2/3}}{3^{2/3}}\ _2F_1 \left (\frac{1}{6}, \frac{2}{3}; \frac{7}{6}; -\frac{1}{3} \right ). \tag6$$

Finally, from DLMF: 15.4.31 we see that $$_2F_1 \left (a, \frac{1}{2} + a; \frac{3}{2}-2a; -\frac{1}{3} \right ) = \left (\frac{8}{9} \right )^{-2a} \frac{\Gamma (\frac{4}{3}) \Gamma (\frac{3}{2} - 2a)}{\Gamma (\frac{3}{2}) \Gamma (\frac{4}{3} - 2a)}.$$ Setting $a = 1/6$ leads to $$_2F_1 \left (\frac{1}{6}, \frac{2}{3}; \frac{7}{6}; -\frac{1}{3} \right ) = \frac{\sqrt[3]{9}}{\sqrt{\pi}} \Gamma \left (\frac{4}{3} \right ) \Gamma \left (\frac{7}{6} \right ).$$ Thus (6) becomes $$_2F_1 \left (1, \frac{2}{3}; \frac{7}{6}; \frac{1}{4} \right ) = \frac{2^{4/3}}{\sqrt{\pi}} \Gamma \left (\frac{4}{3} \right ) \Gamma \left (\frac{7}{6} \right ).$$ On substituting this result into (4), one has \begin{align} I &= \Gamma \left (\frac{2}{3} \right ) \Gamma \left (\frac{4}{3} \right )\\ &= \frac{1}{3} \Gamma \left (\frac{2}{3} \right ) \Gamma \left (\frac{1}{3} \right )\\ &= \frac{1}{3} \Gamma \left (1 - \frac{1}{3} \right ) \Gamma \left (\frac{4}{3} \right )\\ &= \frac{\pi}{3 \sin (\frac{\pi}{3})}\\ &= \frac{2\pi}{3 \sqrt{3}}, \end{align} as required.

The solution heavily exploits symmetry of the integrand.

Let $$I = \int_0^1\frac{x^{2/3}(1-x)^{-1/3}}{1-x+x^2} dx $$ Replace $x$ by $1-x$ and sum up gives $$\tag{1} 2I = \int_0^1 \frac{x^{2/3}(1-x)^{-1/3} + (1-x)^{2/3}x^{-1/3}}{1-x+x^2} dx = \int_0^1 \frac{x^{-1/3}(1-x)^{-1/3}}{1-x+x^2} dx$$

Let $\ln_1$ be complex logarithm with branch cut at positive real axis, while $\ln_2$ be the one whose cut is at negative real axis. Denote $$f(z) = \frac{2}{3}\ln_1(x) - \frac{1}{3}\ln_2 (1-x)$$ Then $f(z)$ is discontinuous along the positive axis, but have different jump in $\arg$ across intervals $[0,1]$ and $[1,\infty)$.

Now integrate $g(z) = e^{f(z)}/(1-z+z^2)$ using keyhole contour. Let $\gamma_1$ be path slightly above $[0,1]$, $\gamma_4$ below. $\gamma_2$ be path slightly above $[1,\infty)$, $\gamma_3$ below. It is easily checked that $$\int_{\gamma 1} g(z) dz = I \qquad \qquad \int_{\gamma 4} g(z) dz = I e^{4\pi i/3}$$ $$\int_{\gamma 2} g(z) dz = e^{\pi i/3} \underbrace{\int_1^\infty \frac{x^{2/3}(x-1)^{-1/3}}{1-x+x^2} dx}_J\qquad \int_{\gamma 3} g(z) dz = e^{\pi i} J$$

If we perform $x\mapsto 1/x$ on $J$, we get $\int_0^1 x^{-1/3}(1-x)^{-1/3}/(1-x+x^2)dx$, thus $J = 2I$ by $(1)$.

Therefore $$I(1-e^{4\pi i/3}) + 2I(e^{\pi i / 3} - e^{\pi i}) = 2\pi i\times \text{Sum of residues of } g(z) \text{ at } e^{\pm 2\pi i /3}$$ From which I believe you can work out the value of $I$.

There is a completely elementary way to solve this. In the end, I do not see how to find the elementary primitive in a simple and intuitive way (if others do, then please edit the answer accordingly), but Rubi helped me. For this reason I post this but make it cw. I also would like to thank @JanG who pointed me to the question and who actually was the one doing the first changes of variables.

Set $$ I=\int_0^1\frac{x^{2/3}(1-x)^{-1/3}}{1-x(1-x)}\,dx. $$ By doing $x\mapsto 1-x$ and adding, others have found that $$ I= \frac{1}{2}\int_0^1\frac{x^{-1/3}(1-x)^{-1/3}}{1-x(1-x)}\,dx. $$ Next, let $x=(1+y)/2$. Then the integral becomes $$ I=2^{2/3}\int_{-1}^1\frac{1}{(1-y^2)^{1/3}(3+y^2)}\,dy= 2^{5/3}\int_{0}^1\frac{1}{(1-y^2)^{1/3}(3+y^2)}\,dy. $$ This is very similar to (and just a $y=\cos t$ away from) the integral @Zacky observes in a comment to the question.

This can be put into Rubi, and surprisinlgy the result is elementary, $$ I=\biggl[\frac{1}{\sqrt{3}}\arctan\Bigl(\frac{\sqrt{3}}{y}\Bigr) +\frac{1}{\sqrt{3}}\arctan\Bigl(\frac{\sqrt{3}\bigl(1-(2-2y^2)^{1/3}\bigr)}{y}\Bigr)-\frac{1}{3}\text{artanh}\,y+\text{artanh}\,\Bigl(\frac{y}{1+(2-2y^2)^{1/3}}\Bigr)\biggr]_0^1 $$ Inserting the boundarys, the upper one gives (using a limit) $2\pi/(3\sqrt{3})$ and the lower one gives $0$. Hence $$ I=\frac{2\pi}{3\sqrt{3}}. $$

Here we piggy back off the solution posted by @pisco, organize the analysis with detail on the definitions of $\arg(z)$ and $\arg(1-z)$, and finish by evaluating the resiudes enclosed by the closed "keyhole contour."

Let $f(z)$ be the function given by

$$f(z)=\frac{z^{2/3}(1-z)^{-1/3}}{z^2-z+1}$$

where choose the branch cut from $0$ to $\infty$ along the positive real axis such that

$$\arg(z)=\begin{cases} 0&, z=x+i0^+\\\\ 2\pi&,z=x+i0^- \end{cases}$$

and we choose the branch cut from $1$ to $\infty$ along the positive real axis with $\arg(1-z)=-\pi+\arg(z-1)$ such that

$$\arg(1-z)=\begin{cases} 0&, 0<x<1\\\\ -\pi&,z=x+i0^+, 1<x\\\\ \pi&, z=x+i0^-, 1<x \end{cases}$$

Then, the integral around the classical "key hole" contour $C$ is

$$\begin{align} \oint_C f(z)\,dz &=(e^{i2(0)/3}e^{-i(0)/3}-e^{i2(2\pi)/3}e^{-i(0)/3})\int_0^1 \frac{x^{2/3}(1-x)^{-1/3}}{x^2-x+1}\,dx\\\\ &+(e^{i2(0)/3}e^{-i(-\pi)/3}-e^{i2(2\pi)/3}e^{-i(\pi)/3})\int_1^\infty \frac{x^{2/3}(x-1)^{-1/3}}{x^2-x+1}\,dx\\\\ &=(1+e^{i\pi/3})\left(\int_0^1 \frac{x^{2/3}(1-x)^{-1/3}}{x^2-x+1}\,dx+\int_1^\infty \frac{x^{2/3}(x-1)^{-1/3}}{x^2-x+1}\,dx\right)\tag1 \end{align}$$

Enforcing the substitution $x\mapsto 1/x$ in the second integral on the right-hand side of $(1)$ reveals

$$\begin{align} \oint_C f(z)\,dz &=(1+e^{i\pi/3})\int_0^1 \frac{x^{2/3}(1-x)^{-1/3}+x^{-1/3}(1-x)^{-1/3}}{x^2-x+1}\,dx\tag2 \end{align}$$

Using the identity $x^{2/3}(1-x)^{-1/3}+x^{-1/3}(1-x)^{2/3}=x^{-1/3}(1-x)^{-1/3}$ and observing that $x^2-x+1=(1-x)^2-(1-x)+1$ we find from $(2)$ that

$$\begin{align} \oint_C f(z)\,dz &=3(1+e^{i\pi/3})\int_0^1 \frac{x^{2/3}(1-x)^{-1/3}}{x^2-x+1}\,dx\\\\ &=3(1+e^{i\pi/3})\int_0^1 \frac{x^{2/3}(1-x)^{-1/3}}{x^2-x+1}\,dx\tag3 \end{align}$$

From the residue theorem we have

$$\begin{align} \oint_C f(z)\,dz&=2\pi i \left(\text{Res}\left(f(z), z=\frac12+i\frac{\sqrt3}2\right)+\text{Res}\left(f(z), z=\frac12-i\frac{\sqrt3}2\right)\right)\\\\ &=2\pi i \left(\frac{e^{i2\pi/9}e^{i\pi/9}}{i2\sqrt 3}+\frac{e^{i10\pi/9}e^{-i\pi/9}}{-i2\sqrt 3}\right)\\\\ &=\frac{2\pi}{\sqrt3} (1+e^{i\pi/3})\tag4 \end{align}$$

Finally, setting $(3)$ and $(4)$ equal yields the coveted result

$$\int_0^1 \frac{x^{2/3}(1-x)^{-1/3}}{x^2-x+1}\,dx=\frac{2\pi }{3\sqrt 3}$$