Let $X\in\{0,1\}^2$ have mean $\mu=\left[\begin{smallmatrix}p_1\\p_2\end{smallmatrix}\right]$ and $\Pr[X_1 = X_2 = 1] = p\le \min\{p_1,p_2\}$. (Note we must have $1-p_1-p_2+p\ge 0$ for the distribution to be well defined.) We can then compute the covariance matrix $\Sigma = E[(X-\mu)(X-\mu)^T] = \left[\begin{smallmatrix}p_1(1-p_1)&p-p_1p_2\\p-p_1p_2&p_2(1-p_2)\end{smallmatrix}\right]$.
I would like to use the Berry Essen bound, and for that we need to upper bound the quantity $\gamma=E[\|\Sigma^{-1/2}(X-\mu)\|_2^3]$.
I believe one should be able to show $$\gamma \le C\left(\tfrac1{\sqrt{p_1(1-p_1)}}+\tfrac1{\sqrt{p_2(1-p_2)}}+\tfrac1{\sqrt{\min\{p_1,p_2\}-p}}\right) , $$ for some universal constant $C>0$.
The symbolic computation of $\Sigma^{-1/2}$ is a bit unwieldy though, and so I wonder if there is some tricks I may use to arrive at this result more neatly?
Or if not, any proof would be appreciated.
Update: Using Pinelli's observation, we can compute $$\begin{align}\gamma |\Sigma|^{3/2} &= p[(1-p_1)(1-p_2)v]^{3/2} \\&+ (p_1-p)[(1-p_1)p_2(1-v)]^{3/2} \\&+ (p_2-p)[p_1(1-p_2)(1-v)]^{3/2} \\&+ (1-p_1-p_2+p)[p_1 p_2 v]^{3/2} , \end{align}$$ where $v = p_1+p_2-2p$ and $|\Sigma|=p_1p_2(1-v)-p^2$.
It seems like I forgot a factor $1/\sqrt{p}$. In fact it seems that bounds of $\gamma \le \frac{2} {\sqrt{p\,(\min\{p_1,p_2\}-p)}} $ or $\gamma \le \frac{1} {\sqrt{p}}+ \frac{2} {\sqrt{\min\{p_1,p_2\}-p}} $ are sufficient and correct.