Integral of squared lower incomplete gamma function

I try to express this integral for the squared lower incomplete gamma function in terms of the same lower incomplete gamma and other common functions $$\int x^{b-1}\gamma(a,x)^2dx\;\;\;(1)$$ When there is no square it is relatively easy to prove $$\int x^{b-1}\gamma(a,x)dx = \frac{1}{b}[x^b\gamma(a,x)-\gamma(a+b,x)]$$ by integrating by parts, starting with change of variables $$u=\gamma(a,x); dv=x^{b-1}dx$$ and then calculating $$\int udv = uv-\int vdu$$ Similar way I'm trying to apply for the first integral (1) by variables change $$u=\gamma(a,x)^2;dv=x^{b-1}dx$$Then I get $$\int x^{b-1}\gamma(a,x)^2dx = \frac{1}{b}[x^b\gamma(a,x)^2-2\int\gamma(a,x)e^{-x}x^{a+b-1}dx]\;\;\;(2)$$Next step again I use change of variables for the integral from right side $$u=\gamma(a,x); dv=e^{-x}x^{a+b-1}dx$$ and then I get$$\int x^{b-1}\gamma(a,x)^2dx = \frac{1}{b}[x^b\gamma(a,x)^2-2(\gamma(a,b)\gamma(a+b,x)-\int\gamma(a+b,x)e^{-x}x^{a-1}dx)]$$ and here I'm stuck. If I try to apply similar integration by parts to the integral from right side I'm getting back the equation (2) Is there a way to get rid of integral on right hand side?

Let’s try a series representation and see if it can be simplified. The series seems to work for a large radius of convergence when integrating by each term:

$$γ(a,z)=-z^a \sum_{n=0}^\infty \frac{(-z)^n}{(a+n)n!},\Gamma(a,z)=\int e^{-x} x^{a-1}dx$$

Therefore: $$\int γ(a+b,x)e^{-x} x^{a-1}dx=-\int e^{-x} x^{a-1}x^{a+b}\sum_{n=0}^\infty \frac{(-1)^nx^n}{(a+b+n)n!}= -\sum_{n=0}^\infty \frac{(-1)^n}{(a+b+n)n!} \int e^{-x}x^{2a+b+n-1}dx= C+\sum_{n=0}^\infty \frac{(-1)^n\Gamma(2a+b+n,x)}{(a+b+n)n!} $$

Here is a technically closed form with the Incomplete Gauss Hypergeometric function which has many papers like these about it, so it is an official special function:

$$\,_2\Gamma_1((a_1,k),a_2;b_1;z)\mathop=^\text{def} \sum_{n=0}^\infty \frac{\Gamma(a_1+n,k)(a_2)_nz^n}{\Gamma(a_1)(b_1)n!}\implies \sum_{n=0}^\infty \frac{\Gamma(a_1+n,k)\Gamma(a_2+n)z^n}{\Gamma(b_1+n)n!}= \frac{\Gamma(a_1)\Gamma(a_2)}{\Gamma(b_1)} \,_2\Gamma_1((a_1,k),a_2;b_1;z),(0)_0\to 1$$

where $(a)_n=\frac{\Gamma(a+n)}{\Gamma(a)}$ is a Pochhammer Symbol:

$$\int γ(a+b,x)e^{-x} x^{a-1}dx= C+\sum_{n=0}^\infty \frac{(-1)^n\Gamma(2a+b+n,x)}{(a+b+n)n!} = C+\sum_{n=0}^\infty \frac{(-1)^n\Gamma(2a+b+n,x) \Gamma(a+b+n)}{\Gamma(a+b+1+n)n!} =\frac{\Gamma(2a+b)}{a+b}\,_2\Gamma_1((2a+b,x),a+b;a+b+1;-1)+C$$ Please correct me and give me feedback!

So, you have reduced your question down to solving the integral $\int\gamma(a+b,x)e^{-x}x^{a-1}dx.$ I'm going to define this integral as $$J(a,b,t) = \int_0^t\gamma(a,x)e^{-x}x^{b}dx,$$ and then your integral is recovered as $J(a+b,b-1,x) + C$.

To solve the integral, I find a recursion relation in $b$ which I solve to find $J$ for all $a$ and $t$, and for integer $b\ge 0$ (note: my $b$ is your $a-1$, apologies for potential confusion here). My method extends straightforwardly to $b\in \mathbb{Z},$ and I'm optimistic that it may be possible to get an analytic extension to all $b$, although I haven't thought too much about this.

You asked in your question, 'when I try to integrate by parts I get the same function back, what can I do'. One way to deal with this is as following. First consider $$\int_0^t\frac{\partial}{\partial x}\left(\gamma(a,x)e^{-x}x^{b}\right)dx = \gamma(a,t)e^{-t}t^{b}.$$ Then expanding the derivative we find the identity that $$\begin{align} \gamma(a,t)e^{-t}t^{b} &= -\int_0^t e^{-2x} x^{a+b-1} dx - \int_0^t \gamma(a,x)e^{-x} x^{b} dx + b \int_0^t \gamma(a,x)e^{-x} x^{b-1} dx\\ &= -\frac{\gamma(a+b,2t)}{2^{a+b}} - J(a,b,t) + b\,J(a,b{-}1,t),\end{align}$$which is roughly what integration by parts tells us about $J$.

Now define $$f(a,b,t) := -\gamma(a,t)e^{-t}t^{b}-\frac{\gamma(a+b,2t)}{2^{a+b}},$$ and we can view this integration by parts identity as a recurrence relation for $J$ in $b$, $$J(a,b,t) = b\,J(a,b{-}1,t) + f(a,b,t).$$ This is a simple linear recurrence relation which can be solved for integer $b\ge 0$ by $$J(a,b,t) = b!\,\left(J(a,0,t) + \sum_{k=0}^b \frac{f(a,k,t)}{k!}\right).$$ Then we need to find $J(a,0,t)$ which is a relatively straightforward integral, I find that $$J(a,0,t) = \gamma(a,t)e^{-t} - \frac{\gamma(a,2t)}{2^a}.$$

Substituting in the functional forms for $J(a,0,t)$ and $f(a,b,t)$ gives a closed form expression for $J(a,b,t)$ for integer $b \ge 0$, which is already more than I expected for this integral. Solving the recursion relation in the opposite direction can also give a similar solution for $b \le 0,$ so this method really solves for $b\in \mathbb{Z}.$ The question then is; can we extend to $b \in \mathbb{R}$ (or $\mathbb{C}$)?

I'm optimistic that this could be possible. First we replace $b!$ with $\Gamma(b+1)$, this is straightforward. We also need to evaluate the sum over $k$ to give a function of $b$ which is valid for $b$ non-integer. There are two terms in the sum and the second one is actually simple to evaluate using the upper incomplete gamma function as an analytic extension. For integer $s$ the incomplete gamma function can be integrated by parts to give $$\Gamma(s,x) = (s-1)!\, e^{-x} \sum_{k=0}^{s-1} \frac{x^k}{k!},$$ so using this identity we can then write $J(a,b,t)$ as $$J(a,b,t) = \Gamma(b+1)\left(\gamma(a,t)e^{-t} - \frac{\gamma(a,2t)}{2^a}\right) - \gamma(a,t)(\Gamma(b+1) - \gamma(b+1,t)) - \frac{b!}{2^a} \sum_{k=0}^{b}\frac{\gamma(a+k,2t)}{k! 2^{k}}. $$

Then to find an expression for your integral for all $b$, we'd need to be able to evaluate $$\sum_{k=0}^{b}\frac{\gamma(a+k,2t)}{k! 2^{k}}$$ to find an expression valid for $b$ non-integer. This kind of analytic extension is exactly the point of defining the gamma function, and it already worked in one of the two sums that came up, so I feel like it could be possible in this case too. I haven't thought about trying to evaluate this sum though, if you needed a closed for expression for all $b$, maybe you could ask another question about evaluating this sum.

(note: I may have made some minus sign errors at different points. They are not crucial to the method.)