We have defined the function $\mathcal{I}:[-1,\infty)^{2}\rightarrow\mathbb{R}$ via the integral representation,
$$\mathcal{I}{\left(a,b\right)}:=\int_{0}^{1}\frac{\ln{\left(1+ax\right)}\ln{\left(1+bx\right)}}{x}\,\mathrm{d}x;~~~\small{\left(a,b\right)\in[-1,\infty)^{2}}.\tag{1}$$
The goal is to find a single closed-form expression for the integral $\mathcal{I}{\left(a,b\right)}$ in terms of elementary functions and polylogarithms that is valid over as large of the domain as possible. The specific definitions for the various polylogarithmic functions employed below have been collected in the appendix at the bottom.
It is obvious from definition $(1)$ that the function $\mathcal{I}{\left(a,b\right)}$ is symmetric under exchange of its arguments, i.e.,
$$\forall \left(a,b\right)\in[-1,\infty)^{2}:\mathcal{I}{\left(a,b\right)}=\mathcal{I}{\left(b,a\right)},\tag{2}$$
so we can assume without loss of generality that $a\le b$.
It is also obvious from definition $(1)$ that the function $\mathcal{I}{\left(a,b\right)}$ will trivially vanish in the special case that one or both of its arguments equals zero:
$$\forall \left(a,b\right)\in[-1,\infty)^{2}:\mathcal{I}{\left(a,0\right)}=\mathcal{I}{\left(0,b\right)}=0.\tag{3}$$
There are additional special cases with simple final expressions. One such special case is that of equal parameters. For $a\in[-1,\infty)$,
$$\begin{align}
\mathcal{I}{\left(a,a\right)}
&=\int_{0}^{1}\frac{\ln^{2}{\left(1+ax\right)}}{x}\,\mathrm{d}x\\
&=2\int_{0}^{1}\frac{\ln^{2}{\left(1+ax\right)}}{2x}\,\mathrm{d}x\\
&=2\,S_{1,2}{\left(-a\right)}.\tag{4}\\
\end{align}$$
Another special case occurs when either of the parameters equals negative unity. For $a\in[-1,\infty)$, we have
$$\begin{align}
\mathcal{I}{\left(a,-1\right)}
&=\int_{0}^{1}\frac{\ln{\left(1+ax\right)}\ln{\left(1-x\right)}}{x}\,\mathrm{d}x\\
&=\int_{0}^{1}\mathrm{d}x\,\frac{\ln{\left(1-x\right)}}{x}\int_{0}^{1}\mathrm{d}y\,\frac{ax}{1+axy}\\
&=\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}y\,\frac{a\ln{\left(1-x\right)}}{1+axy}\\
&=\int_{0}^{1}\mathrm{d}y\int_{0}^{1}\mathrm{d}x\,\frac{a\ln{\left(1-x\right)}}{1+ayx}\\
&=\int_{0}^{1}\mathrm{d}y\int_{0}^{1}\mathrm{d}x\,\frac{1}{1-x}\left[\frac{\ln{\left(1+ayx\right)}}{y}-\frac{\ln{\left(1+ay\right)}}{y}\right];~~~\small{I.B.P.s}\\
&=\int_{0}^{1}\frac{\mathrm{d}y}{y}\int_{0}^{1}\mathrm{d}x\,\frac{\ln{\left(1+ayx\right)}-\ln{\left(1+ay\right)}}{1-x}\\
&=\int_{0}^{1}\frac{\mathrm{d}y}{y}\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(1+ay-ayt\right)}-\ln{\left(1+ay\right)}}{t};~~~\small{\left[1-x=t\right]}\\
&=\int_{0}^{1}\frac{\mathrm{d}y}{y}\int_{0}^{1}\mathrm{d}t\,\frac{\ln{\left(1-\frac{ay}{1+ay}t\right)}}{t}\\
&=\int_{0}^{1}\frac{\mathrm{d}y}{y}\left[-\operatorname{Li}_{2}{\left(\frac{ay}{1+ay}\right)}\right]\\
&=\int_{0}^{1}\frac{\operatorname{Li}_{2}{\left(-ay\right)}+\frac12\ln^{2}{\left(1+ay\right)}}{y}\,\mathrm{d}y;~~~\small{\text{Landen's Identity}}\\
&=\int_{0}^{1}\frac{\operatorname{Li}_{2}{\left(-ay\right)}}{y}\,\mathrm{d}y+\frac12\int_{0}^{1}\frac{\ln^{2}{\left(1+ay\right)}}{y}\,\mathrm{d}y\\
&=\operatorname{Li}_{3}{\left(-a\right)}+S_{1,2}{\left(-a\right)},\tag{5}\\
\end{align}$$
where the so-called dilogarithm identity of Landen mentioned in the above derivation refers to the identity,
$$\operatorname{Li}_{2}{\left(z\right)}+\operatorname{Li}_{2}{\left(\frac{z}{z-1}\right)}+\frac12\ln^{2}{\left(1-z\right)}=0;~~~\small{z\notin[1,\infty)}.\tag{6}$$
Main Result:
Suppose $-1<a<b\land a\neq0\land b\neq0$. Then, using the algebraic identity,
$$2xy=x^{2}+y^{2}-\left(x-y\right)^{2},\tag{7}$$
and the logarithmic identity,
$$\ln{\left(z\right)}-\ln{\left(w\right)}=\ln{\left(\frac{z}{w}\right)};~~~\small{z\ge0\land w\ge0},\tag{8}$$
we have
$$\begin{align}
\mathcal{I}{\left(a,b\right)}
&=\int_{0}^{1}\frac{\ln{\left(1+ax\right)}\ln{\left(1+bx\right)}}{x}\,\mathrm{d}x\\
&=\int_{0}^{1}\frac{2\ln{\left(1+ax\right)}\ln{\left(1+bx\right)}}{2x}\,\mathrm{d}x\\
&=\frac12\int_{0}^{1}\frac{\ln^{2}{\left(1+ax\right)}+\ln^{2}{\left(1+bx\right)}-\left[\ln{\left(1+ax\right)}-\ln{\left(1+bx\right)}\right]^{2}}{x}\,\mathrm{d}x\\
&=\frac12\int_{0}^{1}\frac{\ln^{2}{\left(1+ax\right)}}{x}\,\mathrm{d}x+\frac12\int_{0}^{1}\frac{\ln^{2}{\left(1+bx\right)}}{x}\,\mathrm{d}x\\
&~~~~~-\frac12\int_{0}^{1}\frac{\ln^{2}{\left(\frac{1+ax}{1+bx}\right)}}{x}\,\mathrm{d}x\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-\frac12\int_{0}^{1}\frac{\ln^{2}{\left(\frac{1+ax}{1+bx}\right)}}{x}\,\mathrm{d}x\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-\frac12\int_{1}^{\frac{1+a}{1+b}}\frac{\ln^{2}{\left(y\right)}}{\left(\frac{1-y}{by-a}\right)}\cdot\frac{a-b}{\left(by-a\right)^{2}}\,\mathrm{d}y;~~~\small{\left[\frac{1+ax}{1+bx}=y\right]}\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-\frac12\int_{\frac{1+a}{1+b}}^{1}\frac{\left(b-a\right)\ln^{2}{\left(y\right)}}{\left(1-y\right)\left(by-a\right)}\,\mathrm{d}y\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-\frac12\int_{\frac{1+a}{1+b}}^{1}\frac{\ln^{2}{\left(y\right)}}{1-y}\,\mathrm{d}y-\frac12\int_{\frac{1+a}{1+b}}^{1}\frac{b\ln^{2}{\left(y\right)}}{by-a}\,\mathrm{d}y\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-\frac12\int_{0}^{\frac{b-a}{1+b}}\frac{\ln^{2}{\left(1-t\right)}}{t}\,\mathrm{d}t;~~~\small{\left[1-y=t\right]}\\
&~~~~~-\frac12\int_{\frac{1+a}{1+b}}^{1}\frac{b\ln^{2}{\left(y\right)}}{by-a}\,\mathrm{d}y\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}-\frac12\int_{\frac{1+a}{1+b}}^{1}\frac{b\ln^{2}{\left(y\right)}}{by-a}\,\mathrm{d}y.\tag{9}\\
\end{align}$$
Next, note that $\ln{\left(y\right)}=\ln{\left(\left|y\right|\right)}$ for $y>0$, and $0<\frac{1+a}{1+b}\le y$ for $\left(a,b\right)\in(-1,\infty)^{2}$. Thus,
$$\begin{align}
\mathcal{I}{\left(a,b\right)}
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}-\frac12\int_{\frac{1+a}{1+b}}^{1}\frac{b\ln^{2}{\left(y\right)}}{by-a}\,\mathrm{d}y\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}-\frac12\int_{\frac{1+a}{1+b}}^{1}\frac{b\ln^{2}{\left(\left|y\right|\right)}}{by-a}\,\mathrm{d}y\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}\\
&~~~~~-\frac12\int_{-\frac{\left(1+a\right)b}{a\left(1+b\right)}}^{-\frac{b}{a}}\frac{\ln^{2}{\left(\left|\frac{at}{b}\right|\right)}}{1+t}\,\mathrm{d}t;~~~\small{\left[y=-\frac{at}{b}\right]}\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}\\
&~~~~~-\frac12\int_{-\frac{\left(1+a\right)b}{a\left(1+b\right)}}^{-\frac{b}{a}}\frac{\left[\ln{\left(\left|t\right|\right)}+\ln{\left(\left|\frac{a}{b}\right|\right)}\right]^{2}}{1+t}\,\mathrm{d}t\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}\\
&~~~~~-\frac12\int_{-\frac{\left(1+a\right)b}{a\left(1+b\right)}}^{-\frac{b}{a}}\frac{\ln^{2}{\left(\left|t\right|\right)}+2\ln{\left(\left|\frac{a}{b}\right|\right)}\ln{\left(\left|t\right|\right)}+\ln^{2}{\left(\left|\frac{a}{b}\right|\right)}}{1+t}\,\mathrm{d}t\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}-\frac12\int_{-\frac{\left(1+a\right)b}{a\left(1+b\right)}}^{-\frac{b}{a}}\frac{\ln^{2}{\left(\left|t\right|\right)}}{1+t}\,\mathrm{d}t\\
&~~~~~-\ln{\left(\left|\frac{a}{b}\right|\right)}\int_{-\frac{\left(1+a\right)b}{a\left(1+b\right)}}^{-\frac{b}{a}}\frac{\ln{\left(\left|t\right|\right)}}{1+t}\,\mathrm{d}t-\frac12\ln^{2}{\left(\left|\frac{a}{b}\right|\right)}\int_{-\frac{\left(1+a\right)b}{a\left(1+b\right)}}^{-\frac{b}{a}}\frac{\mathrm{d}t}{1+t}\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}-\frac12\int_{0}^{-\frac{b}{a}}\frac{\ln^{2}{\left(\left|t\right|\right)}}{1+t}\,\mathrm{d}t\\
&~~~~~+\frac12\int_{0}^{-\frac{\left(1+a\right)b}{a\left(1+b\right)}}\frac{\ln^{2}{\left(\left|t\right|\right)}}{1+t}\,\mathrm{d}t-\ln{\left(\left|\frac{a}{b}\right|\right)}\int_{0}^{-\frac{b}{a}}\frac{\ln{\left(\left|t\right|\right)}}{1+t}\,\mathrm{d}t\\
&~~~~~+\ln{\left(\left|\frac{a}{b}\right|\right)}\int_{0}^{-\frac{\left(1+a\right)b}{a\left(1+b\right)}}\frac{\ln{\left(\left|t\right|\right)}}{1+t}\,\mathrm{d}t-\frac12\ln^{2}{\left(\left|\frac{a}{b}\right|\right)}\int_{-\frac{\left(b-a\right)}{a\left(1+b\right)}}^{-\frac{\left(b-a\right)}{a}}\frac{\mathrm{d}u}{u};~~~\small{\left[1+t=u\right]}\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}-\frac12\,\Lambda_{3}{\left(-\frac{b}{a}\right)}+\frac12\,\Lambda_{3}{\left(-\frac{\left(1+a\right)b}{a\left(1+b\right)}\right)}\\
&~~~~~-\ln{\left(\left|\frac{a}{b}\right|\right)}\,\Lambda_{2}{\left(-\frac{b}{a}\right)}+\ln{\left(\left|\frac{a}{b}\right|\right)}\,\Lambda_{2}{\left(-\frac{\left(1+a\right)b}{a\left(1+b\right)}\right)}\\
&~~~~~-\frac12\ln^{2}{\left(\left|\frac{a}{b}\right|\right)}\int_{\frac{1}{1+b}}^{1}\frac{\mathrm{d}v}{v};~~~\small{\left[-\frac{\left(b-a\right)u}{a}=v\right]}\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}-\frac12\left[\Lambda_{3}{\left(-\frac{b}{a}\right)}-\Lambda_{3}{\left(-\frac{\left(1+a\right)b}{a\left(1+b\right)}\right)}\right]\\
&~~~~~-\ln{\left(\left|\frac{a}{b}\right|\right)}\left[\Lambda_{2}{\left(-\frac{b}{a}\right)}-\Lambda_{2}{\left(-\frac{\left(1+a\right)b}{a\left(1+b\right)}\right)}\right]\\
&~~~~~-\frac12\ln^{2}{\left(\left|\frac{a}{b}\right|\right)}\int_{1}^{1+b}\frac{\mathrm{d}w}{w};~~~\small{\left[v=w^{-1}\right]}\\
&=S_{1,2}{\left(-a\right)}+S_{1,2}{\left(-b\right)}-S_{1,2}{\left(\frac{b-a}{1+b}\right)}-\frac12\left[\Lambda_{3}{\left(-\frac{b}{a}\right)}-\Lambda_{3}{\left(-\frac{\left(1+a\right)b}{a\left(1+b\right)}\right)}\right]\\
&~~~~~-\ln{\left(\left|\frac{a}{b}\right|\right)}\left[\Lambda_{2}{\left(-\frac{b}{a}\right)}-\Lambda_{2}{\left(-\frac{\left(1+a\right)b}{a\left(1+b\right)}\right)}\right]-\frac12\ln^{2}{\left(\left|\frac{a}{b}\right|\right)}\ln{\left(1+b\right)}.\tag{10}\\
\end{align}$$
Result $(10)$ provides us with a compact formula for $\mathcal{I}{\left(a,b\right)}$ in terms of logarithms, Kummer polylogarithms, and Nielsen generalized polylogarithms under the parameter conditions, $-1<a<b\land a\neq0\land b\neq0$.
Appendix:
The natural logarithm can be defined for positive real argument via the integral representation,
$$\ln{\left(z\right)}:=\int_{1}^{z}\frac{\mathrm{d}t}{t};~~~\small{z\in\mathbb{R}^{+}}.$$
For positive integer order $n\in\mathbb{N}$, the polylogarithm $\operatorname{Li}_{n}{\left(z\right)}$ may be defined for real argument $z$ recursively via,
$$\begin{cases}
&\operatorname{Li}_{1}{\left(z\right)}=-\ln{\left(1-z\right)};~~~\small{z<1},\\
&\operatorname{Li}_{n}{\left(z\right)}=\int_{0}^{z}\frac{\operatorname{Li}_{n-1}{\left(t\right)}}{t}\,\mathrm{d}t;~~~\small{n\ge2\land z\le1}.\\
\end{cases}$$
An auxiliary function related to the polylogarithm is Kummer's polylogarithm. For positive integer $n\in\mathbb{N}$ and real argument $z\in\mathbb{R}$, we define
$$\Lambda_{n+1}{\left(z\right)}:=\int_{0}^{z}\frac{\ln^{n}{\left(\left|t\right|\right)}}{1+t}\,\mathrm{d}t$$
For positive integer indices $\left(n,p\right)\in\mathbb{N}^{2}$ and real argument $z\in(-\infty,1]$, define the Nielsen generalized polylogarithm via the integral representation,
$$S_{n,p}{\left(z\right)}:=\frac{\left(-1\right)^{n+p-1}}{\left(n-1\right)!\,p!}\int_{0}^{1}\frac{\ln^{n-1}{\left(t\right)}\ln^{p}{\left(1-zt\right)}}{t}\,\mathrm{d}t.$$
In particular,
$$S_{1,p}{\left(z\right)}:=\frac{\left(-1\right)^{p}}{p!}\int_{0}^{1}\frac{\ln^{p}{\left(1-zt\right)}}{t}\,\mathrm{d}t=\frac{\left(-1\right)^{p}}{p!}\int_{0}^{z}\frac{\ln^{p}{\left(1-x\right)}}{x}\,\mathrm{d}x.$$
