It's
$$\int_0^1\frac{\log(1-x+x^2)\log(1+x-x^2)}{x}dx= -2\sum\limits_{k=1}^\infty \frac{(2k-1)!^2}{(4k)!}\sum\limits_{v=0}^{k-1}\frac{1}{k+v}$$
which is $\enspace\approx -0.0848704554500727311… $ .
Already $\enspace\displaystyle -2\sum\limits_{k=1}^{10} \frac{(2k-1)!^2}{(4k)!}\sum\limits_{v=0}^{k-1}\frac{1}{k+v}\enspace$ gives a good approach.
Note: A closed form for such or comparable series is not known to me.
Proof:
$\displaystyle \int_0^1\frac{\log(1-x+x^2)\log(1+x-x^2)}{x}dx=$
$\displaystyle =\int_0^1\lim\limits_{h\to 0}\frac{((1-x+x^2)^h-1)((1+x-x^2)^h-1)}{h^2x}dx$
$\displaystyle =\lim\limits_{h\to 0}\frac{1}{h^2}\int_0^1\frac{((1-x+x^2)^h-1)((1+x-x^2)^h-1)}{x}dx$
$\displaystyle =\lim\limits_{h\to 0}\frac{1}{h^2}\left(\int_0^1\left(\frac{(1-(x-x^2)^2)^h-1}{x}-\frac{(1-x+x^2)^h-1}{x}-\frac{(1-x+x^2)^h-1}{x}\right)dx\right) $
$\displaystyle =\lim\limits_{h\to 0}\frac{1}{h^2}\int_0^1\left(\sum\limits_{k=1}^\infty \binom h k \left(x^{k-1}(-x(1-x)^2)^k -x^{k-1}(-1+x)^k -x^{k-1}(1-x)^k\right) \right) $
$\displaystyle =\lim\limits_{h\to 0}\frac{1}{h^2}\sum\limits_{k=1}^\infty \binom h k \int_0^1\left(x^{k-1}(-x(1-x)^2)^k -x^{k-1}(-1+x)^k -x^{k-1}(1-x)^k\right)dx $
$\displaystyle =\lim\limits_{h\to 0}\frac{1}{h^2}\sum\limits_{k=1}^\infty \binom h k \left((-1)^k\frac{(2k-1)!(2k)!}{(4k)!} -(1+(-1)^k)\frac{(k-1)!k!}{(2k)!}\right) $
$\displaystyle =-\lim\limits_{h\to 0}\frac{1}{h^2}\sum\limits_{k=1}^\infty \left((-1)^{k-1}\binom h k + 2\binom h {2k}\right) \frac{(2k-1)!(2k)!}{(4k)!}$
$\displaystyle =-\sum\limits_{k=1}^\infty \frac{(2k-1)!(2k)!}{(4k)!}\lim\limits_{h\to 0}\frac{1}{h^2}\left((-1)^{k-1}\binom h k + 2\binom h {2k}\right)$
$\displaystyle =-\sum\limits_{k=1}^\infty \frac{(2k-1)!(2k)!}{(4k)!}\frac{H_{2k-1}-H_{k-1}}{k}= -2\sum\limits_{k=1}^\infty \frac{(2k-1)!^2}{(4k)!}\sum\limits_{v=0}^{k-1}\frac{1}{k+v}$
Additional comment:
$$\int_0^1\frac{\log(1-z(x-x^2))\log(1+z(x-x^2))}{x}dx= -2\sum\limits_{k=1}^\infty z^{2k}\frac{(2k-1)!^2}{(4k)!}\sum\limits_{v=0}^{k-1}\frac{1}{k+v}$$
for $\,z\in\mathbb{C}\,$ and $\,|z|\leq 1\,$ .
$$\color{blue}{\int_0^\infty} \color{red}{\int_0^y} f_{T_1, T_2}(x,y) \,\color{red}{dx}\,\color{blue}{dy}$$
$f(x,y)$ is the density function, upon integrating them, we recover the corresponding probability of the interested event.
The integration domain describe the region of interest.
Here $x$ correponds to $T_1$ and $y$ correponds to $T_2$.
The domain of integration does describe the set $\{(x,y)|y \geq x, y\geq 0 , x\geq 0\}$.
First, $y$ can take any nonnegative value. Upon fixing $y$, $x$ can only take value from $0$ up to $y$.
Best Answer
It's integration by parts, i.e. $$\int_0^1 u(x) v'(x) \, dx = [u(x) v(x)]_{x=0}^1 - \int_0^1 u'(x) v(x) \, dx$$ with $u(x) = x$ and $v(x) = (1-x)^{\theta - 1}$.