According to my book (The Nature of computation, page 691):
$$\int_{0}^{2 \pi} \int_{0}^{2 \pi} \ln(\cos^{2} \theta + \cos^{2} \phi) ~d \theta ~d \phi= 16 \pi^2 \left(\frac{C}{\pi}- \frac{\ln2}{2}\right),$$
where $C$ is Catalan's constant. I have tried to derive this expression by looking at integral representations of $C$, but I have not been able to perform the integral. Any help?
Thank you.
Here is one possible reduction that leads to the answer.
Step 1. Let $I$ denote the integral. As in @Takahiro Waki's computation, we utilize several trigonometric identities to write
\begin{align*} I &= \int_{0}^{2\pi}\int_{0}^{2\pi} \log\left( 1 + \frac{\cos2\theta + \cos2\phi}{2} \right) \, \mathrm{d}\theta\mathrm{d}\phi \\ &= \int_{0}^{2\pi}\int_{0}^{2\pi} \log( 1 + \cos(\theta+\phi)\cos(\theta-\phi)) \, \mathrm{d}\theta\mathrm{d}\phi. \end{align*}
Now observe that $(\theta, \phi) \to (\theta-\phi, \theta+\phi) =: (x, y)$, as mapping $(\Bbb{R}/2\pi\Bbb{Z})^2 \to (\Bbb{R}/2\pi\Bbb{Z})^2$ is a 2-1 covering with $\mathrm{d}\theta\mathrm{d}\phi = \frac{1}{2}\mathrm{d}x\mathrm{d}y$. This gives
$$ I = \int_{0}^{2\pi}\int_{0}^{2\pi} \log( 1 + \cos x \cos y ) \, \mathrm{d}x \mathrm{d}y. \tag{1} $$
(See Addendum for a more direct proof.)
Step 2. Applying the McLaurin expansion of the function $z \mapsto \log(1+z)$ to $\text{(1)}$ and integrating term by term, we get
$$ I = -\sum_{n=1}^{\infty} \frac{1}{2n} \left( \int_{0}^{2\pi} \cos^{2n} x \, \mathrm{d}x \right)^2 = -2\pi \sum_{n=1}^{\infty} \frac{\Gamma(n+\frac{1}{2})^2}{n!^2 n}. \tag{2} $$
In order to make a further simplification, notice that for complex $z$ with $|z| < 1$, we have
\begin{align*} \sum_{n=0}^{\infty} (-1)^n \frac{\Gamma(n+\frac{1}{2})}{n!} z^n &= \frac{\sqrt{\pi}}{\sqrt{1+z}}, \tag{3} \\ \sum_{n=1}^{\infty} (-1)^n \frac{\Gamma(n+\frac{1}{2})}{n!n} z^n &= -2\sqrt{\pi} \log\left( \frac{1+\sqrt{1+z}}{2} \right) \tag{4}. \end{align*}
Thus it follows from Parseval's identity that
$$I = 2\pi \int_{-\pi}^{\pi} \log\left( \frac{1+\sqrt{1+e^{i\theta}}}{2} \right) \, \frac{\mathrm{d}\theta}{\sqrt{1 + e^{-i\theta}}}. \tag{5} $$
Now we observe that for $ |\theta| < \pi$, we obtain $\sqrt{1 + e^{-i\theta}} = e^{-i\theta/2}\sqrt{1 + e^{i\theta}}$. Using this, we simplify the above expression as
$$ I = 2\pi \int_{-\pi}^{\pi} \log\left( \frac{1+\sqrt{1+e^{i\theta}}}{2} \right) \frac{e^{i\theta/2}}{\sqrt{1 + e^{i\theta}}} \, \mathrm{d}\theta. $$
Finally, applying the substitution $u = ie^{i\theta/2}$ and replacing the resulting semicircular contour by the linear segment $[-1, 1]$, we find that
$$ I = 4\pi \int_{-1}^{1} \log\left(\frac{1+\sqrt{1-u^2}}{2}\right) \, \frac{\mathrm{d}u}{\sqrt{1-u^2}}. \tag{6} $$
Step 3. It remains to compute $\text{(6)}$. Applying the substitution $u = \sin \theta$, we have
$$ I = 8\pi \int_{0}^{\pi/2} \log\left(\frac{1+\cos\theta}{2}\right) \, \mathrm{d}\theta = 16\pi \int_{0}^{\pi/2} \log\cos(\theta/2) \, \mathrm{d}\theta. $$
The final integral is not hard to compute by using
$$ \log\left|\cos(\theta/2)\right| = \Re\log(1+e^{i\theta}) - \log 2 = -\log 2 + \sum_{n=1}^{\infty} \frac{(-1)^{n-1}}{n} \cos(n\theta), $$
and the result is
$$ I = 16\pi \left( C - \frac{\pi}{2} \log 2 \right) $$
as corrected by @JeanMarie.
Here we collect some arguments which clarifies some steps of the main computation.
Equation (1). Notice that the transform $(x, y) = (\theta-\phi, \theta+\phi)$ maps the square $[0, 2\pi]^2$ to a square $\mathcal{D}$ with vertices $(0, 0)$, $(2\pi, 2\pi)$, $(-2\pi, 2\pi)$ and $(0, 4\pi)$. Now split this square into four non-overlapping parts
$$ \mathcal{D} = \mathcal{D}_1 \cup \mathcal{D}_2 \cup \mathcal{D}_3 \cup \mathcal{D}_4, $$
where
Then by translating each pieces appropriately and reassembling, we find that
Thus utilizing the $2\pi$-periodicity of both $\cos x$ and $\cos y$ we can write
\begin{align*} I &= \frac{1}{2} \iint_{\mathcal{D}} \log( 1 + \cos x \cos y ) \, \mathrm{d}x \mathrm{d}y \\ &= \frac{1}{2} \sum_{i=1}^{4} \iint_{\mathcal{D}_i} \log( 1 + \cos x \cos y ) \, \mathrm{d}x \mathrm{d}y \\ &= 2 \times \frac{1}{2} \iint_{[0, 2\pi]^2} \log( 1 + \cos x \cos y ) \, \mathrm{d}x \mathrm{d}y. \end{align*}
This proves $\text{(1)}$.
Equation (2). This is a simple consequence of the following beta function identity
$$ 2\int_{0}^{\pi/2} \cos^{2s-1}\theta \sin^{2t-1}\theta \, \mathrm{d}\theta = \beta(s, t) = \frac{\Gamma(s)\Gamma(t)}{\Gamma(s+t)}, \quad \Re(s), \Re(t) > 0. $$
Equation (3) and (4). By the generalized binomial theorem, we get
$$ \frac{1}{\sqrt{1+z}} = \sum_{n=0}^{\infty} \binom{-1/2}{n} z^n, \quad |z| < 1. $$
Now by expanding the binomial coefficient, we find that
\begin{align*} \binom{-1/2}{n} &= \frac{(-\frac{1}{2})(-\frac{1}{2}-1)\cdots(-\frac{1}{2}-n+1)}{n!} \\ &= (-1)^n \frac{(n-\frac{1}{2})\cdots(1-\frac{1}{2})}{n!} = (-1)^n \frac{\Gamma(n+\frac{1}{2})}{n!\Gamma(\frac{1}{2})} = (-1)^n \frac{\Gamma(n+\frac{1}{2})}{n!\sqrt{\pi}}. \end{align*}
Plugging this back to the binomial series proves $\text{(3)}$. In order to prove $\text{(4)}$, notice that both sides of $\text{(4)}$ define analytic functions on $|z| < 1$ with value zero at $z = 0$ and that their derivatives coincide:
$$ -2\sqrt{\pi} \frac{\mathrm{d}}{\mathrm{d}z} \log\left( \frac{1+\sqrt{1+z}}{2} \right) = \frac{\sqrt{\pi}}{z}\left( \frac{1}{\sqrt{1+z}} - 1 \right) = \sum_{n=1}^{\infty} (-1)^n \frac{\Gamma(n+\frac{1}{2})}{n!} z^{n-1}. $$
This proves that $\text{(4)}$ is true.
Equation (5). Let $0 < r < 1$. Then using the absolute convergence we can rearrange the sum to write
\begin{align*} &2\pi \log\left( \frac{1+\sqrt{1+re^{i\theta}}}{2} \right)\frac{1}{\sqrt{1+re^{-i\theta}}} \\ &\hspace{9em} = -\sum_{\substack{m \geq 0 \\ n \geq 1}} (-1)^{m+n} \frac{\Gamma(m+\frac{1}{2})\Gamma(n+\frac{1}{2})}{m!n!n} r^{m+n} e^{i\theta(m-n)}. \end{align*}
Now let us integrate both sides with respect to $\theta$ on $[0, 2\pi]$. Since the right-hand side converges uniformly, we can integrate term by term to get
$$ 2\pi \int_{0}^{2\pi} \log\left( \frac{1+\sqrt{1+re^{i\theta}}}{2} \right)\frac{\mathrm{d}\theta}{\sqrt{1+re^{-i\theta}}} = -2\pi \sum_{n=1}^{\infty} \frac{\Gamma(n+\frac{1}{2})^2}{n!^2 n} r^{2n}. $$
As we take limit as $r \uparrow 1$, the left-hand side converges to the left-hand side of $\text{(5)}$ by the dominated convergence theorem. (It is dominated by $C \left| \theta - \pi \right|^{-1/2}$ for some constant $C > 0$.) On the other hand, the right-hand side converges by the monotone convergence theorem to $I$. Therefore $\text{(5)}$ follows.
Let $\cos^2\phi=\frac{(1-r)^2}{4r}$, or, $ \ln r=-2 \sinh^{-1}\cos \phi $, and substitute $x=2\theta$ \begin{align} I=& \int_{0}^{2 \pi} \int_{0}^{2 \pi} \ln(\cos^{2} \theta + \cos^{2} \phi) ~d \theta~ d \phi\\ =& \ 8\int_{0}^{\pi/2} \int_{0}^{\pi} \ln\frac{1+2r\cos x+ r^2}{4r}~dx \ d\phi \end{align} Note that \begin{align} &\int_0^\pi \ln( 1 + 2r\cos x + r^2) dx =-2\sum_{k\ge1}\frac{(-r)^k}k \int_0^\pi \cos kx \ dx =0 \end{align}
$$\int_0^{\frac\pi2} \ln r\ d\phi =-2\int_0^{\frac\pi2}\sinh^{-1}\cos \phi\ d\phi =-2G$$
As a result \begin{align} I=& \ 8\pi\int_{0}^{\pi/2} \ln \frac1{4r}d\phi =16\pi \left(G-\frac\pi2 \ln2\right) \end{align}
$$I=\int_{0}^{2 \pi} \int_{0}^{2 \pi} \ln \left(\frac{\cos2\theta+1}{2} + \frac{\cos2\phi+1}{2}\right) \mathrm d \theta \, \mathrm d \phi.$$
At first $0<\theta+x<2\pi$
$$\left(\int_{0}^{\frac \pi4}\frac{\cos}2 - \int_{\frac \pi4}^{\frac \pi2}\frac{\cos}2-\int_{\frac \pi2}^{\frac {3\pi}4}\frac{\sin}2+\int_{\frac {3\pi}4}^{ {\pi}}\frac{\sin}2 \right) \cdot 2.$$
$$\begin{align} J &=\int_{0}^{2 \pi}\ln \left(\frac{\cos2\theta+1+x}2 \right) \\ &=\int_{0}^{2 \pi}\ln \left(\frac{\sin2\theta+1+x}2 \right) \\ &=\int_{0}^{2 \pi}\ln2+\int_{0}^{2 \pi}\ln(\sin{\theta+1+x})+\int_{0}^{2 \pi}\ln \cos{\theta+1+x} \\ &=4{\pi}\ln2+2J \\ &=-4{\pi}\ln2 \end{align}$$
as same
$$\begin{align} I &=\int_{0}^{2 \pi}\ln \left(\frac{\cos2\theta+1+\phi}2\right) \mathrm d \phi \\ &=-8{\pi}\ln2. \end{align}$$
