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{\vphantom Xx_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_\text{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.