A surprising dilogarithm integral identity arising from a generalised point enclosure problem This question asked:

What is the probability that three points selected uniformly randomly on the unit circle contain a fixed point at distance $x$ from the circle's centre?

I answered that question with a geometric derivation of the following double integral for the desired probability:
$$P(x)=\frac14-\frac3{\pi^3}\int_0^\pi\int_0^a\left(\sin^{-1}\frac{x\sin b}{\sqrt{1+x^2-2x\cos b}}-\sin^{-1}\frac{x\sin a}{\sqrt{1+x^2-2x\cos a}}\right)\,db\,da$$
Below this equation I claimed:

Trying to find a closed-form expression for this is neither simple nor enlightening.

I did give an approximation to $P(x)$ and showed that it had an error of no more than around 0.004 via matplotlib. (In the graphs below, the pole refers to the fixed point.)
$$P(x)\approx\frac{(1-x^2)^{2/3}}4$$

But the relative size of the error irritated me – I had tried more complicated approximations, with only marginal improvements. Then I got the idea of using Taylor series for a better approximation – $P(x)$ is only defined on $[0,1]$, within the radius of convergence for many Taylor series.
So I set to work, using SymPy to compute the Maclaurin (Taylor at $x=0$) series of $P(x)$. After half an hour I got the following:
$$P(x)=\frac14-\frac{3x^2}{2\pi^2}-\frac{3x^4}{8\pi^2}-\frac{x^6}{6\pi^2}-\frac{3x^8}{32\pi^2}-\frac{3x^{10}}{50\pi^2}+\mathcal O(x^{12})$$
I realised that


*

*the initial term and the factors of $\frac3{\pi^2}$ in subsequent terms were echoes of their presence in the expression for $P(x)$, so could be factored out

*computing the $x^{12}$ term would take forever in this form, so I split the integrand into four parts using the logarithmic version of $\sin^{-1}$:
$$\int_0^\pi\int_0^a\left(\sin^{-1}\frac{x\sin b}{\sqrt{1+x^2-2x\cos b}}-\sin^{-1}\frac{x\sin a}{\sqrt{1+x^2-2x\cos a}}\right)\,db\,da=\\
i\left(\int_0^\pi\int_0^a\ln(xe^{ib}-1)\,db\,da-\frac12\int_0^\pi\int_0^a\ln(1+x^2-2x\cos b)\,db\,da\\
-\int_0^\pi\int_0^a\ln(xe^{ia}-1)\,db\,da+\frac12\int_0^\pi\int_0^a\ln(1+x^2-2x\cos a)\,db\,da\right)\tag1$$


With that, I computed the $x^{12}$ term:
$$P(x)=\frac14-\frac3{\pi^2}\left(\frac{x^2}2+\frac{x^4}8+\frac{x^6}{18}+\frac{x^8}{32}+\frac{x^{10}}{50}+\frac{x^{12}}{72}+\mathcal O(x^{14})\right)$$
All the powers were even and had unit fraction coefficients after factoring $\frac3{\pi^2}$ out. I then put the sequence $2,8,18,32,50,72$ into the OEIS and got a surprise… these are twice squares!
Therefore I could write
$$P(x)=\frac14-\frac3{\pi^2}\sum_{n=1}^\infty\frac{x^{2n}}{2n^2}$$
Thus was revealed the series expansion of the dilogarithm, the simplest expression for $P(x)$ and perhaps the most elegant equation I have encountered:
$$\bbox[5px,border:1px solid black]{P(x)=\frac14-\frac3{2\pi^2}\operatorname{Li}_2(x^2)}$$
I confirmed the correctness of this expression by comparing its values to my precomputed samples of $P(x)$, finding nothing but round-off error in the latter:

Of course, such an equality begs explanation, and this is the question I am asking here:

Why is the following identity true?
  $$\int_0^\pi\int_0^a\left(\sin^{-1}\frac{x\sin b}{\sqrt{1+x^2-2x\cos b}}-\sin^{-1}\frac{x\sin a}{\sqrt{1+x^2-2x\cos a}}\right)\,db\,da=\frac\pi2\operatorname{Li}_2(x^2)$$

The logarithmic form $(1)$ of the integrand in $P(x)$ may help. Besides Taylor series, I am open to proofs that exploit the geometry of the situation; here is the diagram I used to justify the integral expression for $P(x)$.

Here $O$ is the origin and the arguments of $A$ and $B$ when treated as complex numbers are $a$ and $b$, lying in $[0,\pi]$ and $[0,a]$ respectively. The argument of the triangle's third point (not shown) lies in $[-\pi,0]$ and has to be between $P$ and $Q$ for the triangle to contain the fixed point $D$ at $(x,0)$. The first term in the integrand $\sin^{-1}\frac{x\sin b}{\sqrt{1+x^2-2x\cos b}}$ is $\angle DOQ$, the second term $\angle DOP$.
 A: Under Construction: This isn't a complete answer, but I wanted to go ahead and point out some symmetries and simplifications to your integral that'll make it much friendlier to work with...

Define the function $f:\left[0,1\right]\times\left[0,\pi\right]\setminus\{\left(1,0\right)\}\rightarrow\left[0,\frac{\pi}{2}\right]$ via the inverse trigonometric expression,
$$f{\left(x,\alpha\right)}:=\arcsin{\left(\frac{x\sin{\left(\alpha\right)}}{\sqrt{1-2x\cos{\left(\alpha\right)}+x^{2}}}\right)},$$
and then define the function $\mathcal{P}:\left[0,1\right]\rightarrow\mathbb{R}$ via the double integral,
$$\begin{align}
\mathcal{P}{\left(x\right)}
&:=\int_{0}^{\pi}\mathrm{d}\alpha\int_{0}^{\alpha}\mathrm{d}\beta\,\left[f{\left(x,\beta\right)}-f{\left(x,\alpha\right)}\right].\\
\end{align}$$

It so happens that any double integral of this form can be reduced to a single-variable integral without even specifying the function $f$ inside the integrand, using basic symmetries.
For any $x\in\left[0,1\right]$, we obtain
$$\begin{align}
\mathcal{P}{\left(x\right)}
&=\int_{0}^{\pi}\mathrm{d}\alpha\int_{0}^{\alpha}\mathrm{d}\beta\,\left[f{\left(x,\beta\right)}-f{\left(x,\alpha\right)}\right]\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\int_{0}^{\alpha}\mathrm{d}\beta\,f{\left(x,\beta\right)}-\int_{0}^{\pi}\mathrm{d}\alpha\int_{0}^{\alpha}\mathrm{d}\beta\,f{\left(x,\alpha\right)}\\
&=\int_{0}^{\pi}\mathrm{d}\beta\int_{0}^{\beta}\mathrm{d}\alpha\,f{\left(x,\alpha\right)}-\int_{0}^{\pi}\mathrm{d}\alpha\int_{0}^{\alpha}\mathrm{d}\beta\,f{\left(x,\alpha\right)}\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\int_{\alpha}^{\pi}\mathrm{d}\beta\,f{\left(x,\alpha\right)}-\int_{0}^{\pi}\mathrm{d}\alpha\int_{0}^{\alpha}\mathrm{d}\beta\,f{\left(x,\alpha\right)}\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left[\int_{\alpha}^{\pi}\mathrm{d}\beta\,f{\left(x,\alpha\right)}-\int_{0}^{\alpha}\mathrm{d}\beta\,f{\left(x,\alpha\right)}\right]\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left[\left(\pi-\alpha\right)f{\left(x,\alpha\right)}-\alpha\,f{\left(x,\alpha\right)}\right]\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left(\pi-2\alpha\right)f{\left(x,\alpha\right)}.\\
\end{align}$$

Given $x\in\left(0,1\right)$, we then find using the trigonometric identity
$$\arcsin{\left(x\right)}=\arctan{\left(\frac{x}{\sqrt{1-x^{2}}}\right)};~~~\small{x^{2}<1},$$
that $\mathcal{P}$ reduces to:
$$\begin{align}
\mathcal{P}{\left(x\right)}
&=\int_{0}^{\pi}\mathrm{d}\alpha\int_{0}^{\alpha}\mathrm{d}\beta\,\left[f{\left(x,\beta\right)}-f{\left(x,\alpha\right)}\right]\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left(\pi-2\alpha\right)f{\left(x,\alpha\right)}\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left(\pi-2\alpha\right)\arcsin{\left(\frac{x\sin{\left(\alpha\right)}}{\sqrt{1-2x\cos{\left(\alpha\right)}+x^{2}}}\right)}\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left(\pi-2\alpha\right)\arctan{\left(\frac{x\sin{\left(\alpha\right)}}{1-x\cos{\left(\alpha\right)}}\right)}\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left(\pi-\alpha\right)\arctan{\left(\frac{x\sin{\left(\alpha\right)}}{1-x\cos{\left(\alpha\right)}}\right)}\\
&~~~~~-\int_{0}^{\pi}\mathrm{d}\alpha\,\alpha\arctan{\left(\frac{x\sin{\left(\alpha\right)}}{1-x\cos{\left(\alpha\right)}}\right)}\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left(\pi-\alpha\right)\arctan{\left(\frac{x\sin{\left(\alpha\right)}}{1-x\cos{\left(\alpha\right)}}\right)}\\
&~~~~~-\int_{0}^{\pi}\mathrm{d}\alpha\,\left(\pi-\alpha\right)\arctan{\left(\frac{x\sin{\left(\alpha\right)}}{1+x\cos{\left(\alpha\right)}}\right)};~~~\small{\left[\alpha\mapsto\pi-\alpha\right]}\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left(\pi-\alpha\right)\left[\arctan{\left(\frac{x\sin{\left(\alpha\right)}}{1-x\cos{\left(\alpha\right)}}\right)}-\arctan{\left(\frac{x\sin{\left(\alpha\right)}}{1+x\cos{\left(\alpha\right)}}\right)}\right]\\
&=\int_{0}^{\pi}\mathrm{d}\alpha\,\left(\pi-\alpha\right)\arctan{\left(\frac{x^{2}\sin{\left(2\alpha\right)}}{1-x^{2}\cos{\left(2\alpha\right)}}\right)}\\
&=\frac14\int_{0}^{2\pi}\mathrm{d}\varphi\,\left(2\pi-\varphi\right)\arctan{\left(\frac{x^{2}\sin{\left(\varphi\right)}}{1-x^{2}\cos{\left(\varphi\right)}}\right)};~~~\small{\left[2\alpha=\varphi\right]}\\
&=\frac14\int_{0}^{2\pi}\mathrm{d}\varphi\,\left(-\varphi\right)\arctan{\left(\frac{x^{2}\sin{\left(\varphi\right)}}{1-x^{2}\cos{\left(\varphi\right)}}\right)};~~~\small{\left[\varphi\mapsto2\pi-\varphi\right]}\\
&=\frac14\int_{0}^{2\pi}\mathrm{d}\varphi\,\left(\pi-\varphi\right)\arctan{\left(\frac{x^{2}\sin{\left(\varphi\right)}}{1-x^{2}\cos{\left(\varphi\right)}}\right)}\\
&=\frac12\int_{0}^{\pi}\mathrm{d}\varphi\,\left(\pi-\varphi\right)\arctan{\left(\frac{x^{2}\sin{\left(\varphi\right)}}{1-x^{2}\cos{\left(\varphi\right)}}\right)}\\
&=\frac12\int_{0}^{\pi}\mathrm{d}\varphi\int_{\varphi}^{\pi}\mathrm{d}\omega\,\arctan{\left(\frac{x^{2}\sin{\left(\varphi\right)}}{1-x^{2}\cos{\left(\varphi\right)}}\right)}\\
&=\frac12\int_{0}^{\pi}\mathrm{d}\omega\int_{0}^{\omega}\mathrm{d}\varphi\,\arctan{\left(\frac{x^{2}\sin{\left(\varphi\right)}}{1-x^{2}\cos{\left(\varphi\right)}}\right)}...\\
\end{align}$$

...TBC
A: Continuing from the last expression for the original double integral in David's answer:
$$\mathcal P(x)=\int_0^\pi\int_0^a(\dots)\,db\,da=\frac12\int_0^\pi\int_0^\omega\tan^{-1}\frac{x^2\sin\varphi}{1-x^2\cos\varphi}\,d\varphi\,d\omega$$
Using the logarithmic form of $\tan^{-1}$ we can rewrite the integrand of the RHS as
$$\begin{align}\tan^{-1}\frac{x^2\sin\varphi}{1-x^2\cos\varphi}&=\frac i2\ln\frac{1-\frac{ix^2\sin\varphi}{1-x^2\cos\varphi}}{1+\frac{ix^2\sin\varphi}{1-x^2\cos\varphi}}\\
&=\frac i2\ln\frac{\frac{1-x^2\cos\varphi-ix^2\sin\varphi}{1-x^2\cos\varphi}}{\frac{1-x^2\cos\varphi+ix^2\sin\varphi}{1-x^2\cos\varphi}}\\
&=\frac i2\ln\frac{1-x^2(\cos\varphi+i\sin\varphi)}{1-x^2(\cos\varphi-i\sin\varphi)}\\
&=\frac i2\ln\frac{1-x^2e^{i\varphi}}{1-x^2e^{-i\varphi}}\\
&=\frac i2(\ln(1-x^2e^{i\varphi})-\ln(1-x^2e^{-i\varphi}))
\end{align}$$
Using the Maclaurin series of $\ln(1-z)$, which is valid since $|x^2e^{i\varphi}|=x^2\le1$:
$$\begin{align}
\frac i2(\ln(1-x^2e^{i\varphi})-\ln(1-x^2e^{-i\varphi}))&=\frac i2\left(\sum_{n=1}^\infty-\frac{(x^2e^{i\varphi})^n}n-\sum_{n=1}^\infty-\frac{(x^2e^{-i\varphi})^n}n\right)\\
&=-\frac i2\sum_{n=1}^\infty\frac{x^{2n}}n(e^{in\varphi}-e^{-in\varphi})\\
&=-\frac i2\sum_{n=1}^\infty\frac{x^{2n}}n\cdot2i\sin n\varphi\\
&=\sum_{n=1}^\infty\frac{x^{2n}}n\sin n\varphi
\end{align}$$
Therefore
$$\frac12\int_0^\pi\int_0^\omega\tan^{-1}\frac{x^2\sin\varphi}{1-x^2\cos\varphi}\,d\varphi\,d\omega=\frac12\int_0^\pi\int_0^\omega\sum_{n=1}^\infty\frac{x^{2n}}n\sin n\varphi\,d\varphi\,d\omega$$
We can move the summation and $\frac{x^{2n}}n$ factor to the outside, then perform simple termwise integration to get the desired result:
$$\begin{align}
\frac12\int_0^\pi\int_0^\omega\sum_{n=1}^\infty\frac{x^{2n}}n\sin n\varphi\,d\varphi\,d\omega&=\frac12\sum_{n=1}^\infty\frac{x^{2n}}n\int_0^\pi\int_0^\omega\sin n\varphi\,d\varphi\,d\omega\\
&=\frac12\sum_{n=1}^\infty\frac{x^{2n}}n\int_0^\pi\left[-\frac1n\cos n\varphi\right]_0^\omega\,d\omega\\
&=\frac12\sum_{n=1}^\infty\frac{x^{2n}}n\int_0^\pi\frac1n(1-\cos n\omega)\,d\omega\\
&=\frac12\sum_{n=1}^\infty\frac{x^{2n}}{n^2}\int_0^\pi(1-\cos n\omega)\,d\omega\\
&=\frac12\sum_{n=1}^\infty\frac{x^{2n}}{n^2}\left[\omega-\frac1n\sin n\omega\right]_0^\pi\\
&=\frac12\sum_{n=1}^\infty\frac{x^{2n}}{n^2}\cdot\pi\\
&=\frac\pi2\sum_{n=1}^\infty\frac{x^{2n}}{n^2}\\
\mathcal P(x)&=\frac\pi2\operatorname{Li}_2(x^2)
\end{align}$$

More generally, we have the following identity for $|z|\le1$:
$$\int_0^\pi\int_0^\omega\tan^{-1}\frac{z\sin\varphi}{1-z\cos\varphi}\,d\varphi\,d\omega=\pi\operatorname{Li}_2(z)$$
