Let \begin{align} G(x_1,x_2,x_3)=\int_0^{2\pi}\int_0^{2\pi}\int_0^{2\pi} \frac{d\theta_1 d\theta_2 d\theta_3}{(2\pi)^3} p(\vec{x},\vec{\theta}) \log\left( p(\vec{x}, \vec{\theta})\right) \end{align} be defined on the probability simplex $\sum_{i=1}^{3} x_i =1$ and $x_i \geq 0$, where \begin{align} p(\vec{x},\vec{\theta}) = 1 + 2 \sum_{i=1}^3 \sum_{j=1}^i \sqrt{x_i x_j} \cos(\theta_i- \theta_j) \end{align} Note $p(\vec{x},\vec{\theta}) \geq 0$ (point wise i.e. for any $\vec{x}, \vec{\theta}$ in the domain) and it can be shown that $G \geq 0$.

I'd like to understand where $G$'s maximum is and if there is a simple closed-form expression for the maximum value.

It can be seen numerically that $x_1 = x_2 = x_3 = 1/3$ is the unique maximum point which yields $G_{max} \approx 0.329524$. My idea is to try to show $G$ is strictly concave, which I also see to be numerically true. But this seems harder than it looks... More generally, any help/advice on any method to proceed analytically would be much appreciated, thanks.