Double integral of symmetric polylogarithmic function over rectangular region

This question was inspired by M.N.C.E.'s wonderful response here. While exploring the possibility of generalizing his result, I found that a significant part of the problem reduced to evaluating the following integral, which is interesting in its own right.

Problem: Define $\mathcal{I}{\left(a,b\right)}$ for $a,b\in\mathbb{R}^{+}\land 0<a,b\le1$ by the double integral $$\mathcal{I}{\left(a,b\right)}:=\iint_{[0,a]\times[0,b]}\mathrm{d}x\mathrm{d}y\,\frac{\chi_{2}{\left(x\right)}-\chi_{2}{\left(y\right)}}{x-y}.\tag{1}$$ Find a closed form evaluation of $\mathcal{I}{\left(a,b\right)}$ in terms of polylogarithms.

Note that $\chi_{\nu}{\left(z\right)}$ represents the Legendre chi function of order $\nu$, and may be defined via the series representation

$$\chi_{\nu}{\left(z\right)}:=\sum_{k=0}^{\infty}\frac{z^{2k+1}}{\left(2k+1\right)^{\nu}};~~~\small{\left[\left|z\right|<1\right]}.\tag{2}$$

Alternatively, the Legendre chi function may be expressed in terms of polylogarithms as

$$\chi_{\nu}{\left(z\right)}=\frac{\operatorname{Li}_{\nu}{\left(z\right)}-\operatorname{Li}_{\nu}{\left(-z\right)}}{2}.\tag{3}$$

(Disclaimer: I think I have nearly solved the above problem, and I plan to post and accept my own answer if I can finish the last bits still troubling me. That said, if anyone else happens to post a solution before I do, I will happily accept that over my own solution. :) )