Preliminary. When $a,b>\frac12$, $c\not=0, -1,-2, ...$, one have (using Mellin transform)
- $I=I(a,b,c)=\int_0^{\infty } \, _2F_1(a,b;c;-x){}^2 \, dx=\frac{\Gamma (c)^2 G_{4,4}^{3,3}\left(1\left| \begin{array}{c} 0,1-a,1-b,c-1 \\ 0,a-1,b-1,1-c \\ \end{array} \right.\right)}{\Gamma (a)^2 \Gamma (b)^2}$
When $a+b<c+\frac12$ it equals
- $I=\frac{(c-1) \, _4F_3(1,a,b,2-c;2-a,2-b,c;1)}{(a-1) (b-1)}-\frac{2 \pi ^3 \csc (2 \pi a) \csc (2 \pi b) \Gamma (c)^2 \cos (\pi (a+b)) \Gamma (a+b-1) \Gamma \left(-a-b+c+\frac{1}{2}\right)}{\Gamma \left(\frac{3}{2}-a\right) \Gamma (a)^2 \Gamma \left(\frac{3}{2}-b\right) \Gamma (b)^2 \Gamma \left(c-\frac{1}{2}\right) \Gamma (c-a) \Gamma (c-b)}$.
Examples: Special cases of formula above. When $a+b<\frac32$:
- $\int_0^{\infty } \, _2F_1(a,b;1;-x){}^2 \, dx=-\frac{\pi ^{3/2} 2^{-2 a-2 b+3} \csc (2 \pi a) \csc (2 \pi b) \cos (\pi (a+b)) \Gamma \left(-a-b+\frac{3}{2}\right) \Gamma (a+b-1)}{\Gamma (2-2 a) \Gamma (a)^2 \Gamma (2-2 b) \Gamma (b)^2}$
When $c>3/2$:
- $\int_0^{\infty } \, _2F_1(a,2-a;c;-x){}^2 \, dx=-\frac{\frac{\sec (\pi a) \Gamma (c)^2}{\Gamma (c-a) \Gamma (a+c-2)}+(c-1)^2}{(a-1)^2 (2 c-3)}$
By using analytic continuation more results are found. Here are $2$ more examples:
- $\int_0^{\infty } \, _2F_1\left(\frac{11}{8},\frac{5}{8};\frac{5}{4};-x\right){}^2 \, dx=\frac{8}{9}-\frac{256 \sqrt[4]{2} \Gamma \left(\frac{9}{8}\right) \Gamma \left(\frac{5}{4}\right)}{9 \sqrt{\left(2-\sqrt{2}\right) \pi } \Gamma \left(-\frac{1}{8}\right)}$
- $\int_0^{\infty } \, _2F_1\left(\frac{3}{4},\frac{5}{6};1;-x\right){}^2 \, dx=-\frac{2 \sqrt[6]{2} \left(\sqrt{3}-3\right) \pi \Gamma \left(\frac{7}{12}\right)}{\Gamma \left(\frac{5}{12}\right) \Gamma \left(\frac{2}{3}\right) \Gamma \left(\frac{3}{4}\right)^2}$
Question: What more can be found for this kind of integrals? This is rather an open question, and any suggestions will be appreciated.
If you see this as a Mellin transform with $s=1$, then by the Ramanujan master theorem you are talking about the integral being $$ \int_0^\infty x^{s-1} \;_2F_1(a,b;c;-x)^2 \; dx = \Gamma(s) C_{-s} $$ for the power series parameterised as $$ _2F_1(a,b;c;-x)^2 = \sum_{k=0}^\infty \frac{(-1)^k}{k!}C_k x^k $$ but it depends on the RMT still holding for this product of power series. I find this is why the negative argument in the hypergeometric function works well from $(-1)^k x^k = (-x)^k$. So perhaps think about the Cauchy product $$ \left(\sum_{i=0}^\infty \frac{(a)_i (b)_i}{(c)_i i!} (-x)^i\right)\left(\sum_{j=0}^\infty \frac{(a)_j (b)_j}{(c)_j j!} (-x)^j\right) $$ also you might want to rewrite the $\pi\csc(\pi s)$ terms in the form $\Gamma(s)\Gamma(1-s)$ to spot patterns. I have a note about spotting patterns in the gamma functions, I'll see if it can find it...
Edit: The following might be useful if there is a way to consider a 'confluence' of sorts from integrals of the form $$ \int_0^\infty \int_0^\infty x_1^{s_1-1} x_2 ^{s_2-1} f_1(x_1) f_2(x_2) \; dx_1 dx_2 \to \int_0^\infty x^{s-1} f(x) f(x) \; dx $$ we could view this as a multidimensional Mellin transform, but I have found there can be different results from the order of integration. If some Fubini type condition is there, then:
If functions $f_k(x)$ have Mellin transforms $g_k(s)$ and coefficients are included then this nested $D$ dimensional type Mellin transform of the product of the functions is given by $$ \mathcal{M}_D\left[\prod_{k=1}^n f_k\left(\alpha_k \prod_{l=1}^n x_l^{a_{kl}}\right) \right] = \frac{\prod_{k=1}^n \alpha_k^{-(A^\top)^{-1}_k \mathbf{s}}}{|\det(A)|}\prod_{k=1}^n g_k((A^\top)^{-1}_k \mathbf{s}) $$ where $A_{kl}=a_{kl}$.
An example Solve $$ I = \int_0^\infty \int_0^\infty \int_0^\infty x_1^{s_1-1} x_2^{s_2-1} x_3^{s_3-1} e^{-\frac{\alpha x_1 x_2}{x_3}}J_n(\beta x_1^2 x_2)\mathrm{Ai}(\gamma x_3) \; dx_1 dx_2 dx_3 $$ with Bessel function $J_n(x)$, Airy function $\mathrm{Ai}(x)$. We have that $f_1(x) = e^{-x}$,$f_2(x) = J_n(x)$, $f_3(x) = \mathrm{Ai}(x)$. We look up that $$ g_1(s) = \Gamma(s) $$ $$ g_2(s) = \frac{2^{s-1} \Gamma \left(\frac{n}{2}+\frac{s}{2}\right)}{\Gamma \left(\frac{n}{2}-\frac{s}{2}+1\right)} $$ $$ g_3(s) = \frac{3^{\frac{2 s}{3}-\frac{7}{6}} \Gamma \left(\frac{s}{3}+\frac{1}{3}\right) \Gamma \left(\frac{s}{3}\right)}{2 \pi } $$ we inspect the integrand and find the coefficient matrix $$ A = \begin{bmatrix} 1 & 1 & -1 \\ 2 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \;\; (A^\top)^{-1} = \begin{bmatrix} -1 & 2 & 0 \\ 1 & -1 & 0 \\ -1 & 2 & 1 \end{bmatrix}, \;\; \det(A) = -1 $$ we have $$ I = \alpha^{s_1 - 2 s_2}\beta^{s_2-s_1}\gamma^{s_1-2 s_2-s_3} \Gamma(2s_2-s_1) \frac{2^{s_1-s_2-1} \Gamma \left(\frac{n}{2}+\frac{s_1-s_2}{2}\right)}{\Gamma \left(\frac{n}{2}-\frac{s_1-s_2}{2}+1\right)} \frac{3^{\frac{2 (2s_2-s_1+s_3)}{3}-\frac{7}{6}} \Gamma \left(\frac{(2s_2-s_1+s_3)}{3}+\frac{1}{3}\right) \Gamma \left(\frac{(2s_2-s_1+s_3)}{3}\right)}{2 \pi } $$
Implication When I see your example result, there are immediate patterns that hint at this linear combination of variables type approach $$ A = -\frac{\pi ^{3/2} 2^{-2 a-2 b+3} \csc (2 \pi a) \csc (2 \pi b) \cos (\pi (a+b)) \Gamma \left(-a-b+\frac{3}{2}\right) \Gamma (a+b-1)}{\Gamma (2-2 a) \Gamma (a)^2 \Gamma (2-2 b) \Gamma (b)^2} $$ for example $2^{-2 a - 2 b + 3}=4^{-a-b+3/2}$ and the linear combination is $-a-b+3/2$ as seen in the gamma function. One possible goal is to split your expression apart into the product of $N$ distinct Mellin transforms and reverse engineer the original integral as a product of simpler integrals?
We can rewrite your first result using $$ \cos\left(\frac{\pi s}{2}\right) = \frac{\pi}{\Gamma\left(\frac{1}{2} + \frac{s}{2}\right)\Gamma\left(\frac{1}{2}-\frac{s}{2}\right)} $$ and $$ \pi \csc(\pi s) = \Gamma(s)\Gamma(1-s) $$ to $$ A = -\pi\csc (2 \pi a) \pi\csc (2 \pi b) \cos (\frac{\pi}{2} (2a+2b)) \frac{4^{-a-b+3/2}}{\pi^{1/2}}\frac{\Gamma \left(-a-b+\frac{3}{2}\right) \Gamma (a+b-1)}{\Gamma (2-2 a) \Gamma (a)^2 \Gamma (2-2 b) \Gamma (b)^2} $$ $$ A = - \pi^{1/2} 4^{-a-b+3/2}\frac{\Gamma(2a)\Gamma(1-2a) \Gamma(2b)\Gamma(1-2b)}{\Gamma\left(\frac{1}{2} + a+b\right)\Gamma\left(\frac{1}{2}-a-b\right)} \frac{\Gamma \left(-a-b+\frac{3}{2}\right) \Gamma (a+b-1)}{\Gamma (2-2 a) \Gamma (a)^2 \Gamma (2-2 b) \Gamma (b)^2} $$
As a base observation $$ \int_0^\infty \int_0^\infty x_1^{s_1-1} x_2 ^{s_2-1} \;_2F_1(a,b;c;-x_1)\;_2F_1(a,b;c;-x_2) \; dx_1 dx_2 = \frac{\Gamma (c)^2 \Gamma (\text{s1}) \Gamma (\text{s2}) \Gamma (a-\text{s1}) \Gamma (a-\text{s2}) \Gamma (b-\text{s1}) \Gamma (b-\text{s2})}{\Gamma (a)^2 \Gamma (b)^2 \Gamma (c-\text{s1}) \Gamma (c-\text{s2})} $$ I feel there might be a parameterisation with $s_1 =a+b-1$ and $s_2=3/2-a-b$, therefore $\Gamma(1/2+a+b) = \Gamma(3/2+s_1)$ and $\Gamma(1/2-a-b)=\Gamma(s_2-1)$. But I have fried my brain for the time being...