We may use the following integral representation for the incomplete gamma function (DLMF ref. see identity 8.6.4):
$$\Gamma{\left(\alpha,z\right)}=\frac{z^{\alpha}e^{-z}}{\Gamma{\left(1-\alpha\right)}}\int_{0}^{\infty}\frac{t^{-\alpha}e^{-t}}{z+t}\,\mathrm{d}t;~~~\small{\left|\arg{\left(z\right)}\right|<\pi,~\Re{\left(\alpha\right)}<1}.$$
Given $n\in\mathbb{N}^{+}\land a,b,\epsilon\in\mathbb{R}\land1<a\land1\le b\le n\land0<\epsilon$, define $I_{n}{\left(a,b;\epsilon\right)}$ via the integral
$$I_{n}{\left(a,b;\epsilon\right)}=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{2ie^{-ibx}\left(bx-i\right)\left[a\left(-ix\right)^{a}\,\Gamma{\left(-a,-ix\right)}\right]^{n}}{x^{2}+\epsilon^{2}}.$$
Then,
$$\begin{align}
I_{n}{\left(a,b;\epsilon\right)}
&=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{2ie^{-ibx}\left(bx-i\right)\left[a\left(-ix\right)^{a}\,\Gamma{\left(-a,-ix\right)}\right]^{n}}{x^{2}+\epsilon^{2}}\\
&=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{2ie^{-ibx}\left(bx-i\right)}{x^{2}+\epsilon^{2}}\left[\frac{e^{ix}}{\Gamma{\left(a\right)}}\int_{0}^{\infty}\frac{t^{a}e^{-t}}{t-ix}\,\mathrm{d}t\right]^{n}\\
&=\frac{1}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{2ie^{-ibx}\left(bx-i\right)}{x^{2}+\epsilon^{2}}e^{inx}\left[\int_{0}^{\infty}\frac{t^{a}e^{-t}}{t-ix}\,\mathrm{d}t\right]^{n}\\
&=\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(bx-i\right)e^{i\left(n-b\right)x}}{x^{2}+\epsilon^{2}}\left[\int_{0}^{\infty}\frac{t^{a}e^{-t}}{t-ix}\,\mathrm{d}t\right]^{n}\\
&=\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(bx-i\right)e^{i\left(n-b\right)x}}{x^{2}+\epsilon^{2}}\prod_{k=1}^{n}\int_{0}^{\infty}\mathrm{d}t_{k}\,\frac{t_{k}^{a}e^{-t_{k}}}{t_{k}-ix}\\
&=\small{\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(bx-i\right)e^{i\left(n-b\right)x}}{x^{2}+\epsilon^{2}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,e^{-\sum_{k=1}^{n}t_{k}}\prod_{k=1}^{n}\left(\frac{t_{k}^{a}}{t_{k}-ix}\right)}\\
&=\small{\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\left(\prod_{k=1}^{n}t_{k}^{a}\right)e^{-\sum_{k=1}^{n}t_{k}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(bx-i\right)e^{i\left(n-b\right)x}}{\left(x^{2}+\epsilon^{2}\right)\prod_{k=1}^{n}\left(t_{k}-ix\right)}}\\
&=\small{\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\frac{\left(\prod_{k=1}^{n}t_{k}^{a}\right)}{\exp{\left(\sum_{k=1}^{n}t_{k}\right)}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(-bx-i\right)e^{-i\left(n-b\right)x}}{\left(x^{2}+\epsilon^{2}\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}}\\
&=:\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\frac{\left(\prod_{k=1}^{n}t_{k}^{a}\right)}{\exp{\left(\sum_{k=1}^{n}t_{k}\right)}}\,f_{n}{\left(b;\epsilon\right)},\\
\end{align}$$
where for $n\in\mathbb{N}^{+}\land b,\epsilon\in\mathbb{R}\land1\le b\le n\land0<\epsilon$ we've defined the auxiliary function denoting the innermost integration,
$$f_{n}{\left(b;\epsilon\right)}:=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(-bx-i\right)e^{-i\left(n-b\right)x}}{\left(x^{2}+\epsilon^{2}\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}.$$
For the integration over $x$, I will appeal to a well-suited proposition from Gradshteyn's Tables. Now, I usually try to make my posts as self-contained as possible by proving any non-trivial lemmas I plan on using, rather than taking the lazy way out and simply citing the result from an outside source. But I can't resist doing so, partly out of laziness, but mostly just to show off how Gradshteyn can on occasion be eerily clairvoyant.
Gradshteyn 3.386: Given the conditions
$$-1<\Re{\left(\nu_{0}\right)}\land0<\Re{\left(\beta_{k}\right)}\land\sum_{k=0}^{n}\Re{\left(\nu_{k}\right)}<1\land0<p,$$
we have the following two results:
$$\int_{-\infty}^{\infty}\frac{\left(ix\right)^{\nu_{0}}e^{-ipx}\prod_{k=1}^{n}\left(\beta_{k}+ix\right)^{\nu_{k}}}{\beta_{0}-ix}\mathrm{d}x=2\pi e^{-\beta_{0}p}\beta_{0}^{\nu_{0}}\prod_{k=1}^{n}\left(\beta_{0}+\beta_{k}\right)^{\nu_{k}},$$
and
$$\int_{-\infty}^{\infty}\frac{\left(ix\right)^{\nu_{0}}e^{-ipx}\prod_{k=1}^{n}\left(\beta_{k}+ix\right)^{\nu_{k}}}{\beta_{0}+ix}\mathrm{d}x=0.$$
The following partial fraction decomposition is easily verified:
$$\frac{-bx-i}{x^{2}+\epsilon^{2}}=\frac{1}{2i\epsilon}\left[\frac{\left(1-b\epsilon\right)}{\epsilon-ix}+\frac{\left(1+b\epsilon\right)}{\epsilon+ix}\right].$$
Then, assuming $0<\Re{\left(t_{k}\right)}\land b<n$,
$$\begin{align}
f_{n}{\left(b;\epsilon\right)}
&=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(-bx-i\right)e^{-i\left(n-b\right)x}}{\left(x^{2}+\epsilon^{2}\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}\\
&=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{1}{2i\epsilon}\left[\frac{\left(1-b\epsilon\right)}{\epsilon-ix}+\frac{\left(1+b\epsilon\right)}{\epsilon+ix}\right]\frac{e^{-i\left(n-b\right)x}}{\prod_{k=1}^{n}\left(t_{k}+ix\right)}\\
&=\frac{1-b\epsilon}{2i\epsilon}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{e^{-i\left(n-b\right)x}}{\left(\epsilon-ix\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}\\
&~~~~~+\frac{1+b\epsilon}{2i\epsilon}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{e^{-i\left(n-b\right)x}}{\left(\epsilon+ix\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}\\
&=\frac{1-b\epsilon}{2i\epsilon}\cdot\frac{2\pi e^{-\left(n-b\right)\epsilon}}{\prod_{k=1}^{n}\left(\epsilon+t_{k}\right)}.\\
\end{align}$$
Thus,
$$\begin{align}
I_{n}{\left(a,b;\epsilon\right)}
&=\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\frac{\left(\prod_{k=1}^{n}t_{k}^{a}\right)}{\exp{\left(\sum_{k=1}^{n}t_{k}\right)}}\,f_{n}{\left(b;\epsilon\right)}\\
&=\frac{2\pi\left(1-b\epsilon\right)e^{-\left(n-b\right)\epsilon}}{\left[\Gamma{\left(a\right)}\right]^{n}\epsilon}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\frac{\left(\prod_{k=1}^{n}t_{k}^{a}\right)}{\exp{\left(\sum_{k=1}^{n}t_{k}\right)}}\cdot\frac{1}{\prod_{k=1}^{n}\left(\epsilon+t_{k}\right)}\\
&=\frac{2\pi\left(1-b\epsilon\right)e^{-\left(n-b\right)\epsilon}}{\left[\Gamma{\left(a\right)}\right]^{n}\epsilon}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\prod_{k=1}^{n}\left(\frac{t_{k}^{a}e^{-t_{k}}}{\epsilon+t_{k}}\right)\\
&=\frac{2\pi\left(1-b\epsilon\right)e^{-\left(n-b\right)\epsilon}}{\left[\Gamma{\left(a\right)}\right]^{n}\epsilon}\left[\int_{0}^{\infty}\mathrm{d}t\,\left(\frac{t^{a}e^{-t}}{\epsilon+t}\right)\right]^{n}\\
&=\frac{2\pi\left(1-b\epsilon\right)e^{-\left(n-b\right)\epsilon}}{\left[\Gamma{\left(a\right)}\right]^{n}\epsilon}\left[\epsilon^{a}e^{\epsilon}\,\Gamma{\left(1+a\right)}\,\Gamma{\left(-a,\epsilon\right)}\right]^{n}\\
&=2\pi a^{n}\left(1-b\epsilon\right)\epsilon^{na-1}e^{b\epsilon}\left[\Gamma{\left(-a,\epsilon\right)}\right]^{n}.\\
\end{align}$$
Interesting question. We may start from the definition of the Beta function:
$$B(m, n) = \int_0^1 t^{m-1}(1-t)^{n-1}\ \text{d}t$$
and rewrite it when $m\to 2m$:
$$B(2m, n) = \int_0^1 t^{2m-1}(1-t)^{n-1}\ \text{d}t$$
$$B(2m, n) = \int_0^1 t^{2m-1} \frac{(1-t)^n}{1-t}\ \text{d}t$$
Now, since the range of integration is $[0, 1]$, we are allowed to make use of the geometric series
$$\frac{1}{1-t} = \sum_{k = 0}^{+\infty} t^k$$
Hence
$$B(2m, n) = \int_0^1 t^{2m-1} (1-t)^n \sum_{k = 0}^{+\infty} t^k\ \text{d}t = \sum_{k = 0}^{+\infty} \int_0^1 t^{2m-1} (1-t)^n t^k\ \text{d}t$$
And easily write:
$$\sum_{k = 0}^{+\infty} \int_0^1 t^{2m-1+k} (1-t)^n\ \text{d}t$$
Calling now
$$2m-1+k = a ~~~~~~~~~~~ n = b$$
We notice that the integral is well known:
$$ \int_0^1 t^a (1-t)^b\ \text{d}t \equiv B(a+1, b+1)$$
Then we end up with the partial result (re-expanding $a$ and $b$):
$$B(2m, n) = \sum_{k = 0}^{+\infty} B(2m+k, n+1)$$
That series does exist and it does converge to a known result:
$$\sum_{k = 0}^{+\infty} B(2m+k, n+1) = \frac{\Gamma (n) \Gamma (2 m+n+1) B(2 m,n+1)}{\Gamma (n+1) \Gamma (2 m+n)}$$
What you end up with is a sort of recursive relation for the Beta function:
$$B(2m, n) = \frac{\Gamma (n) \Gamma (2 m+n+1) B(2 m,n+1)}{\Gamma (n+1) \Gamma (2 m+n)}$$
BUT
The above expression can be simplified!
$$\frac{\Gamma (n) \Gamma (2 m+n+1) B(2 m,n+1)}{\Gamma (n+1) \Gamma (2 m+n)} \equiv \frac{\Gamma (2 m) \Gamma (n)}{\Gamma (2 m+n)}$$
What we obtained is actually nothing than what we would have obtained by simply substituting at the beginning $m\to 2m$ in the Gamma function / Beta function definition.
$$B(2m, n) = \frac{\Gamma (2 m) \Gamma (n)}{\Gamma (2 m+n)}$$
This really suggest that such a particular duplication formula for the Beta function may not exist at all, since all you need is the Gamma function and ITS duplication formula, through which you can evaluate $\Gamma(2m)$.
Seems like that this is the only "duplication formula" for the beta function.
(Also, I found nothing on reviews or literature).
Best Answer
Using the hypergeometric representation of the incomplete Beta function \begin{equation} B_x\left( c,d \right)=\frac{x^c}{c}{}_2F_1\left( c,1-d;1+c;x \right) \end{equation} the integral can be written as \begin{align} I\left( a,b,c,d \right)&=\int_{0}^1x^{a-1}(1-x)^{b-1}B_x(c,d)\,dx\\ &=\frac{1}{c}\int_{0}^1x^{a+c-1}(1-x)^{b-1}{}_2F_1\left( c,1-d;1+c;x \right)\,dx\\ &=\frac{1}{c}\sum_{k=0}^\infty \frac{(c)_k(1-d)_k}{(1+c)_kk!} \int_0^1x^{a+c+k-1}(1-x)^{b-1}\,dx\\ &=\frac{1}{c}\sum_{k=0}^\infty \frac{(c)_k(1-d)_k}{(1+c)_kk!} \frac{\Gamma(b)\Gamma(a+c+k)}{\Gamma(a+b+c+k)}\\ &=\frac{1}{c}\frac{\Gamma(b)\Gamma(a+c)}{\Gamma(a+b+c)}\sum_{k=0}^\infty \frac{(c)_k(1-d)_k}{(1+c)_kk!} \frac{(a+c)_k}{(a+b+c)_k}\\ &=\frac{1}{c}\frac{\Gamma(b)\Gamma(a+c)}{\Gamma(a+b+c)}\,{}_3F_2\left( 1-d,a+c,c;1+c,a+b+c ;1\right) \end{align} Using this identity for the generalized hypergeometric function: \begin{equation} {}_3F_2\left( a_1,a_2,a_3;b_1,b_2;1 \right)=\frac{\Gamma(b_1)\Gamma(b_1+b_2-a_1-a_2-a_3)}{\Gamma(b_1-a_1)\Gamma(b_1+b_2-a_2-a_3)}{}_3F_2\left( a_1,b_2-a_2,b_2-a_3;b_2,b_1+b_2-a_2-a_3;1 \right) \end{equation} Here we choose $a_1=1-d,a_2=a+c,a_3=c,b_1=a+b+c,b_2=1+c$ to obtain \begin{align} I\left( a,b,c,d \right)&=\frac{1}{c}\frac{\Gamma(b)\Gamma(a+c)}{\Gamma(a+b+c)}\frac{\Gamma(a+b+c)\Gamma(b+d)}{\Gamma(a+b+c+d-1)\Gamma(b+1)}\,{}_3F_2\left( 1-d,1-a,1;1+c,1+b ;1\right)\\ &=\frac{1}{bc}\frac{\Gamma(a+c)\Gamma(b+d)}{\Gamma(a+b+c+d-1)}\,{}_3F_2\left( 1,1-a,1-d;1+b,1+c ;1\right)\\ &=\frac{a+b+c+d-1}{bc}B(a+c,b+d)\,{}_3F_2\left( 1,1-a,1-d;1+b,1+c ;1\right) \end{align} which gives simple results if $a$ or $d$ are positive integer, as the hypergeometric series is finite: as $1-a\le0$ or $1-b\le0$, in the series definition of the hypergeometric function, there are $\operatorname{min}(a,b)$ terms (numerator of the coefficients cancel after that). One can also check the given result when $a=c$ and $b=d$.