Let $\boldsymbol{X} \sim Mult_K(n, \boldsymbol{p})$ and $k, l \in \{1, \dots, K\}$.
I need $Cov(X^2_k, X^2_l)$, but I am not sure how to calculate $E(X^2_k \cdot X^2_l)$. With the law of iterated expectation I got
$E[X^2_{k} \cdot X^2_{l}] = E[E[X^2_{k}X^2_{l}|X_{k}]] = E[X^2_{k} \cdot E[X^2_{l}|X_{k}]]$.
Then I use $X_{l}|X_{k} \sim binomial(n-X_k, p_l)$ (not sure if this is correct).
That gives me $E[X^2_{l}|X_{k}] = p_l^2X^2_{k} + [-p_l(1-p_l+2np_l)]X_{k} + np_l(1-p_l) +n^2p^2_l$ and I substitute this into $E[X^2_{k} \cdot E[X^2_{l}|X_{k}]]$.
After few steps I got
\begin{align*} E[X^2_k \cdot X^2_l] &= - \sigma_{kl} \bigg[-1 + n + p_l \cdot (n^2 - 3n + 2) + p_k \cdot (n^2 - 4n + 3) \\ &\qquad \qquad + p_l p_k \cdot[-6n^2 + 4n - 10+ n(n-1)^2] \\ &\qquad \qquad - p^2_k \cdot 2(n^2-3n+2) + p^2_k p_l \cdot 2(-2n^3+13n^2-25n+14) \\ &\qquad\qquad + p_k^3 p_l \cdot (n^3-6n^2+11n-5)\bigg] \end{align*}
I would expect that the result will be symetric for $k$ and $l$ and the high power of $p_k$ seems to be wrong too. I would rather express the $E[X^2_{l}|X_{k}]$ with the correlation coefficient but I do not know how. So I will appreciate any suggestion - where I am wrong or how to do this easier.