Methods to attack integrals that include $(1+x)^{a}\ln^{b}(1+x)$ in the integrand
I am looking for systematic methods to attack the following class of integrals involving logarithmic functions
$$\begin{aligned} I_{0} &= \int_{0}^{1}(1+x)^{a}\ln^{m}(1+x)\,\mathrm{d}x \\ I_{1} &= \int_{0}^{1} x^{a}(1+x)^{b}\ln^{m}x\ln^{n}(1+x)\,\mathrm{d}x \\ I_{2} &= \int_{0}^{1} (1-x)^{a}(1+x)^{b}\ln^{m}(1-x)\ln^{n}(1+x)\,\mathrm{d}x \\ I_{3} &= \int_{0}^{1}x^{a}(1-x)^{b}(1+x)^{c}\ln^{\ell}x\ln^{m}(1-x)\ln^{n}(1+x)\,\mathrm{d}x\end{aligned}$$
where $a, b, c$ are real numbers (usually integers and half-integers) and $\ell, m, n$ are whole numbers such that the integrals converge. Of course, the first three are special cases of $I_{3}$, but I thought it would be best to tackle the easier cases first before moving to the fully general case in case different techniques are needed.
Motivation. The beta function can be written in the form
$$ \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} = \int_{0}^{1}t^{a-1}(1-t)^{b-1}\,\mathrm{d}t. $$
From here, it is straightforward to compute integrals of the type
$$ \int_{0}^{1}x^{a}(1-x)^{b}\ln^{m}x\ln^{n}(1-x)\,\mathrm{d}x $$
through various series expansions or derivatives of the beta function. Thus one is not afraid of any powers on the logarithms, as one merely has to keep enough terms in the expansion, or compute enough derivatives.
Attempt at $I_{0}$. The presence of $1+x$ as a power and in the logarithm means that one cannot use the beta function directly. I am able to get away with it with the integral $I_{0}$ for now by shifting the boundaries and proceeding as
$$\begin{aligned} \int_{0}^{1}(1+x)^{a+\epsilon}\,\mathrm{d}x &= \sum_{m=0}^{\infty}\frac{\epsilon^{m}}{m!}\int_{0}^{1}(1+x)^{a}\ln^{m}(1+x)\,\mathrm{d}x \\ &= \int_{1}^{2}x^{a+\epsilon}\,\mathrm{d}x = \frac{2^{a+1+\epsilon} - 1}{a+1+\epsilon} = \frac{2^{a+1}e^{\epsilon\ln 2} - 1}{a+1+\epsilon} \end{aligned}$$
and keeping the appropriate order coefficient by expanding the exponential in its Taylor series and the denominator in a power series. For example,
$$\begin{aligned} \int_{0}^{1}(1+x)^{1+\epsilon}\,\mathrm{d}x &= \sum_{m=0}^{\infty}\frac{\epsilon^{m}}{m!}\int_{0}^{1}(1+x)\ln^{m}(1+x)\,\mathrm{d}x = \frac{4e^{\epsilon\ln 2} - 1}{2+\epsilon} \\ &\approx \left[4\left(1 + \epsilon\ln 2 + \frac{1}{2}\ln^{2}2\,\epsilon^{2}\right) - 1\right]\frac{1}{2}\left(1 - \frac{\epsilon}{2} + \frac{\epsilon^{2}}{4}\right) \\ &\approx \frac{1}{2}\left(3 + 4\ln 2\,\epsilon + 2\ln^{2}2\,\epsilon^{2}\right)\left(1 - \frac{\epsilon}{2} + \frac{\epsilon^{2}}{4}\right) \end{aligned}$$
so if one requires the $\epsilon^{2}$ coefficient, he merely keep all the terms, multiply by $2! = 2$ to account for the factorial in the original expansion, and arrive at
$$ \int_{0}^{1}(1+x)\ln^{2}(1+x)\,\mathrm{d}x = \frac{3}{4} - 2\ln 2 + 2\ln^{2}2 $$
with the integrals
$$ \int_{0}^{1}(1+x)\ln(1+x)\,\mathrm{d}x = -\frac{3}{4} + 2\ln 2 \quad\text{and}\quad \int_{0}^{1}(1+x)\,\mathrm{d}x = \frac{3}{2}$$
found for "free" if we were not lazy with the multiplication. So we have a method for doing $I_{0}$ for any powers desired.
Attempts at $I_{1}$. The above technique does not seem to generalize to the other integrals. For example, I am currently trying to derive the result
$$ \int_{0}^{1}x(1+x)\ln^{2}x\ln(1+x)\,\mathrm{d}x = -\frac{1}{4}\zeta(3) + \frac{4}{27}\ln 2 - \frac{5}{36}\zeta(2) + \frac{37}{72} \tag{1}$$
and many other integrals of a similar form.
(Edit: a method for $(1)$ has been found, see $(2)$ below)
$$ \int_{0}^{1}(1+x)\ln x\ln^{2}(1+x)\,\mathrm{d}x = -\frac{1}{8}\zeta(3) + 5\ln 2 - 2\ln^{2}2 + \frac{1}{4}\zeta(2) - \frac{23}{8} \tag{2}$$
The class of integrals represented by $I_{1}$ seem difficult, as the zeta functions are familiar from the beta function, the $\ln 2$ is familiar from the treatment of $I_{0}$, but how an answer like this combines these constants is elusive. Hypergeometric functions seem out of reach, as I seem to lose information about the answer from the Pochhammer symbols. To demonstrate this, one may derive from the Euler integral representation of ${}_{2}F_{1}$
$$\begin{aligned} \int_{0}^{1}x^{\epsilon}(1+x)^{\delta}\,\mathrm{d}x &= \frac{1}{1+\epsilon}\,{}_{2}F_{1}(-\delta, 1+\epsilon; 2+\epsilon; -1) \\ \int_{0}^{1}(1-x)^{\epsilon}(1+x)^{\delta}\,\mathrm{d}x &= \frac{1}{1+\epsilon}\,{}_{2}F_{1}(-\delta, 1; 2+\epsilon; -1) \\ \int_{0}^{1}x^{\epsilon}(1-x)^{\delta}(1+x)^{\lambda}\,\mathrm{d}x &= \frac{\Gamma(1+\epsilon)\Gamma(1+\delta)}{\Gamma(2+\delta+\epsilon)}\,{}_{2}F_{1}(-\lambda, 1+\epsilon; 2+\delta+\epsilon; -1). \end{aligned}$$
To take the first one as an example, we can write this expression as
$$ \frac{1}{1+\epsilon}\,{}_{2}F_{1}(-\delta, 1+\epsilon; 2+\epsilon; -1) = \sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!}\frac{(-\delta)_{k}}{1+k+\epsilon} = \sum_{k=0}^{\infty}{\delta\choose k}\frac{1}{1+k+\epsilon}.$$
Thus one can use power series again to find the desired coefficient of $\epsilon$, but the $(-\delta)_{k}$ makes things difficult. If we wanted to find the coefficient of $\epsilon\delta^{2}$, for example as needed in $(2)$, we would need to sum up an infinite series, since every term $k \geq 2$ will contribute to this coefficient. This particular series does numerically converge to the correct answer, but I cannot find a way to relate it to well-known constants such as $\zeta(3)$ and $\ln 2$.
I have also searched up some other integrals of similar forms on this site, and have found that many of the answers proceed with pages of horrendous algebra involving double series, harmonic numbers, etc. which I would like to avoid if possible as I am not very familiar with manipulating harmonic numbers as they appear in series. Of course, I will accept manipulations of harmonic numbers, etc. if it is the only straightforward way to proceed.
Lastly, the references below may be of some use. The second reference may well be the key to solving this problem but I would like some clarification as dealing with the $\ln(1+x)$ still seems problematic...
Devoto, A. and Duke, D. Table of integrals and formulae for Feynman diagram calculations. La Rivista del Nuovo Cimento, 1984.
Vermaseren, J. A. M. Harmonic sums, Mellin transforms and integrals. International Journal of Modern Physics A, 2012.
Sofo, A. Integrals of logarithmic and hypergeometric functions. Communications in Mathematics, 2016.
Solution 1:
Here are some thoughts about $I_2$. The main objective is to "reduce" the product of the logartihms to a single one in a place, I will try to show how this works by an example. Since I saw this in your answer:
$$ I=\int_{0}^{1}\ln^{2}(1-x)\ln(1+x)\,\mathrm{d}x = \frac{7}{2}\zeta(3) - 6 + 2\zeta(2) + \frac{2}{3}\ln^{3}2 - 2\ln^{2}2 - 2\zeta(2)\ln 2 + 4\ln 2$$ I'll try to start with it. If we substract $(a-b)^3$ from $(a+b)^3$ we will remain only with the terms that has an odd power in $(a-b)^3$ since the rest of them will get canceled by the ones from $(a+b)^3$. This is better visualised looking at the Pascal's triangle. $$(a+b)^3-(a-b)^3=2(3a^2b+b^3)\Rightarrow a^2b=\frac16\left((a+b)^3-(a-b)^3-2b^3\right)$$ Choosing $a=\ln(1-x)$ and $b=\ln(1+x)$ we get: $$I=\frac16 \int_0^1\left(\ln^3(1-x^2)-\ln^3\left(\frac{1-x}{1+x}\right)-2\ln^3(1+x)\right)$$ Now by letting $x^2=t$ in the first integral, $\frac{1-x}{1+x}=t$ we get: $$I=\frac1{12}\int_0^1 t^{-\frac12}\ln^3(1-t)dt +\frac13\int_0^1 \frac{\ln^3t}{(1+t)^2}dt-\frac13 \int_0^1 \ln^3(1+t)dt$$ Now to get the result we can see that the first one is doable by beta function, on the second one we can expand the denominator into power series and the third one is elementary.
As another example, consider: $$I=\int_0^1 \frac{\ln^4 (1+x) +6\ln^2(1+x)\ln^2(1-x)}{x}dx$$ This time we can use another heuristic relation: $$(a+b)^4+(a+b)^4=2(a^4+6a^2b^2+b^4)\Rightarrow a^4+6a^2b^2=\frac12((a+b)^4 +(a-b)^4)-b^4$$ $$I=\frac12 \int_0^1 \frac{\ln^4 (1-x^2)}{x}dx +\frac12 \int_0^1 \frac{\ln^4\left(\frac{1-x}{1+x}\right)}{x}dx-\int_0^1 \frac{\ln^4(1-x)}{x}dx$$ In the first integral let $x^2= t$ to obtain: $$\int_0^1 \frac{\ln^4 (1-x^2)}{x}dx=\frac12 \int_0^1 \frac{\ln^4(1-x)}{x}dx =\frac12 \int_0^1 \frac{\ln^4 x}{1-x}dx$$ For the second integral substitute as usual: ${\frac{1-x}{1+x}=t}$ $$\int_0^1 \frac{\ln^4\left(\frac{1-x}{1+x}\right)}{x}dx=2\int_0^1 \frac{\ln^4 t}{1-t^2}dt$$ And for the last one we just do a $1-x\rightarrow x$ and combine with the first one, since they are identical. $$\Rightarrow I=\int_0^1 \frac{\ln^4 x}{1-x^2}dx -\frac34 \int_0^1 \frac{\ln^4 x}{1-x}dx$$ $$I= \sum_{n=0}^\infty \int_0^1 x^{2n} \ln^4 x dx -\frac34 \sum_{n=0}^\infty \int_0^1 x^n \ln^4 x dx$$ $$I=24\sum_{n=0}^\infty \frac{1}{(2n+1)^5}-18 \sum_{n=0}^\infty \frac{1}{(n+1)^5}=\frac{21}{4}\zeta(5)$$ Of course this example was set up to work-out nice and we don't always have this luxury.
But as seen a start for $I_2$ is first to try and rewrite $a^m a^n$ into a form of: $(a-b)^v\pm (a+b)^w +\dots$ where $a=\ln(1-x),\, b=\ln(1+x)$.
It appears that this method works for $I_3$ too, an example is found here, and the idea remains the same: Keep $\ln^ lx$ in a place and rewrite $\ln^m(1-x) \ln^n(1+x)$ as a sum or difference.
Here is an interesting integral that is linked to many other integrals such as this one that also contain other methods to attack those integrals in your question.
Solution 2:
For $$J_0=\int(1+x)^{a}\log^{m}(1+x)\,dx$$ let $1+x=e^t$ to end with $$J_0=\int e^{(a+1) t} t^m\,dt=-t^{m+1} E_{-m}(-(a+1) t)$$ where appears the exponential integral function.
So, $$I_0=\int_0^1(1+x)^{a}\log^{m}(1+x)\,dx=\int_1^{\log(2)}e^{(a+1) t} t^m\,dt$$ and using the relations between the exponential integral function and the complete/incomplete gamma functions $$I_0=(-1)^{m+1}\frac{ \Gamma (m+1)-\Gamma (m+1,-(a+1) \log (2))}{(a+1)^{m+1} }$$
Solution 3:
This is a partial answer; I am still looking for methods to evaluate $I_{1}$ that includes $\ln^{n}(1+x)$ with a power $n \geq 2$, and $I_{2}, I_{3}$ in general.
I have been able to "grasp some straws" and evaluate some integrals in $I_{1}$, specifically of the form
$$ \int_{0}^{1}x^{a}(1+x)^{b}\ln^{m}x\ln(1+x)\,\mathrm{d}x$$
where we can see that the power on the $\ln(1+x)$ is $1$. This implies that it is possible to derive the result $(1)$ as stated in the question. I shall use $(1)$ as an example to illustrate the technique.
We begin by considering the integral
$$\begin{aligned} \int_{0}^{1}x^{1+\epsilon}(1+x)^{1+\delta}\,\mathrm{d}x &= \frac{1}{2+\epsilon}\,{}_{2}F_{1}(-1-\delta,2+\epsilon;3+\epsilon;-1) \\ &= \sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!}\frac{(-1-\delta)_{k}}{2+k+\epsilon}.\end{aligned}$$
It turns out that it is essential to separate the $k=0$ and $k=1$ contributions from the rest of the sum. The reason is because when we manipulate $(-1-\delta)_{k}/k!$, we will find that $(k-1)(k)$ in the denominator does not cancel with the terms in $(-1-\delta)_{k}$. It is easy to see that the $(1+x)^{b}$ term is responsible for this. Therefore we write
$$\begin{aligned}\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!}\frac{(-1-\delta)_{k}}{2+k+\epsilon} = \frac{1}{2+\epsilon} + \frac{1+\delta}{3+\epsilon} + \sum_{k=2}^{\infty}\frac{(-1)^{k}}{2+k+\epsilon}\frac{(-1-\delta)_{k}}{(1)_{k}}. \end{aligned}$$
Since we are interested in finding integrals with the power on $\ln(1+x)$ being $1$, we are interested in the $\delta$ coefficient, specifically $\epsilon^{2}\delta$. Therefore the first term cannot contribute to the integral. The second term is trivial by expanding into a power series, and we obtain a contribution of $1/27$ from it. The rest of the sum is as follows:
$$\begin{aligned} \sum_{k=2}^{\infty}\frac{(-1)^{k}}{2+k+\epsilon}\frac{(-1-\delta)_{k}}{(1)_{k}} &= \sum_{k=2}^{\infty}\frac{(-1)^{k}}{2+k+\epsilon}\frac{(-1-\delta)(-\delta)(1-\delta)\cdots(-2-\delta+k)}{(1)(2)\cdots(k-2)(k-1)(k)} \\ &= \sum_{k=2}^{\infty}\left[\frac{(-1)^{k}}{2+k+\epsilon}\frac{\delta(\delta+1)}{k(k-1)}(1-\delta)\left(1-\frac{\delta}{2}\right)\cdots\left(1-\frac{\delta}{k-2}\right)\right]\end{aligned}$$
Then the $\epsilon^{2}\delta$ coefficient can be obtained in a similar way as what we did for $I_{0}$. The $\epsilon^{2}\delta$ coefficient of the sum is then
$$ \sum_{k=2}^{\infty}\frac{(-1)^{k}}{k(k-1)}\sum_{j=0}^{\infty}\frac{(-1)^{j}\epsilon^{j}}{(2+k)^{j+1}} \to \sum_{k=2}^{\infty}\frac{(-1)^{k}}{k(k-1)}\frac{1}{(2+k)^{3}}.$$
The series is easily done by partial fractions, as each of the decomposed series converges separately. Taking care of signs, the shifting of indices, the $1/27$ out front, the overall factor of $2!1! = 2$, and remembering $\eta(s) = (1-2^{1-s})\zeta(s)$, $\eta(1) = \ln 2$, the result in $(1)$ is found.
This technique relied on the fact that a factor of $\delta$ could be pulled out in front of the sum so that the product $\prod_{n=1}^{k-2}\left(1-\delta/n\right)$ could not contribute to the integral. Because of this, I do not see how this technique can generalize for $\ln^{n}(1+x)$, $n \geq 2$. Nor do I see a way to compute $I_{2}$ or $I_{3}$ because the Pochhammer symbols with $\epsilon$ do not cancel. So I wonder even more how something like
$$ \int_{0}^{1}\ln^{2}(1-x)\ln(1+x)\,\mathrm{d}x = \frac{7}{2}\zeta(3) - 6 + 2\zeta(2) + \frac{2}{3}\ln^{3}2 - 2\ln^{2}2 - 2\zeta(2)\ln 2 + 4\ln 2$$
is found.