$$
\begin{align}
(x-\mu)^T\Sigma^{-1}(x-\mu) &= x^T\Sigma^{-1}x + \mu^T\Sigma^{-1}\mu - \mu^T\Sigma^{-1}x - x^T\Sigma^{-1}\mu \\
&= x^T\Sigma^{-1}x + \mu^T\Sigma^{-1}\mu - 2 x^T\Sigma^{-1}\mu \\
&=x\Sigma^{-1/2}\Sigma^{-1/2}x + \mu\Sigma^{-1/2}\Sigma^{-1/2}\mu - 2 \left<x\Sigma^{-1/2}|\Sigma^{-1/2}\mu\right> \\
&=:||a||^2+||b||^2-2\left<a|b\right>\text{, where } a=\Sigma^{-1/2}x, b=\Sigma^{-1/2}\mu
\end{align}$$
In the second step we used the fact that $\mu^T\Sigma^{-1}x$ is a real number, so it is equal to its transpose, and the fact that $\Sigma$ has to be symmetric.
In a sense, the result is of the form "two squared terms minus twice an inner product", similar to the Euclidean distance form.
Also I think your Euclidean distance formula has a mistake (or at least we're using different notations). It should be $I_{xy}=\sqrt{||x||^2+||y||^2-2\left<x|y\right>}$.
I noticed the mistake in the Euclidean distance formula also, but it comes from a paper, and I thought it was a misunderstanding from myself.
– user12910 Oct 10 '22 at 08:47