I know how I want to answer this question. I'll write up the easy parts here, and leave the hard part for you :).
First some minor changes. It will be convenient to clear out denominators and work with $\log \left( 4 + e^{2 \pi i x} + e^{- 2 \pi i x} + e^{2 \pi i y} + e^{- 2 \pi i y} \right)$. That just changes the constant term of your Fourier series by $\log 2$. Next, it is convenient to focus on $$ \int_0^1 \int_0^1 \log \left( 4 + e^{2 \pi i x} + e^{- 2 \pi i x} + e^{2 \pi i y} + e^{- 2 \pi i y} \right) e^{2 \pi i m x} e^{2 \pi i n y} dx dy.$$ A simple linear transformation goes between this and the cosine formulation. Let $S = \{ (z,w) : |z|=|w|=1 \}$. So we are dealing with $$\frac{1}{(2 \pi i )^2} \int_S \log \left( 4+z+z^{-1} + w +w^{-1} \right) z^{m-1} w^{n-1} dz dw.$$ Dropping out the $4 \pi ^2$, we want to show the integrand is of the form $a \pi + b$.
UPDATE: Thanks to fedja for pointing out that I had oversimplified the next paragraph.
Assuming that $(m,n) \neq (0,0)$, we can integrate by parts with respect to one of the two variables, let's say $z$. Once we do that, we will have a quantity of the form $$ (\mbox{rational number}) \cdot \int_S \frac{(z-z^{-1}) w^k z^{\ell} dw dz} {4+w+w^{-1}+z+z^{-1}}$$
So we'd like to show this quantity is of the form $a+b \pi$.
As fedja points out, we need to be careful here. Without the $z-z^{-1}$ term, the integral diverges like $\int \int dt du/(t^2 + u^2)$$\int \int ds dt/(s^2 + t^2)$ near $(-1,-1)$.
Whew! Now comes the actual hard part. Let
$$E:=\{ (z,w) \in (\mathbb{C}^*)^2 : \ 4+z+z^{-1}+w+w^{-1} =0 \}.$$
ThisThis is an elliptic curve with four punctures. As Bjorn points out, this is an elliptic curve with four puncturesa nodal cubic and can be parameterized as
$$(z,w) = \left( \frac{1-u}{u(1+u)}, \frac{u(u-1)}{1+u} \right).$$
We'll come back to this point later.
The $2$-form $dw dz/(4+z+z^{-1}+w+w^{-1})$ has a simple pole on $E$. Let $\omega$ be the $1$-form on $E$ which is the residue of that $2$-form.
I think there should be a curve $\gamma$ in $E$ such that $S$ is homotopic, in $(\mathbb{C}^*)^2 \setminus E$, to a tubular neighborhood of $\gamma$. So $$\int_S \frac{w^k z^{\ell} (z-z^{-1}) dw dz} {4+w+w^{-1}+z+z^{-1}} = \int_{\gamma} \omega w^k z^{\ell} (z - z^{-1}) .$$
And, well, this is where I wimp out. If $E$ is any elliptic curve defined over $\mathbb{Q}$we substitute in the above parameterization, $\gamma$this will be the integral around a curveclosed loop of some rational function in $E$$\mathbb{Q}(u)$. In particular, we can compute this integral by residues and $\eta$ any meromorphic $1$-form defined over $\mathbb{Q}$, then $\int_{\gamma} \eta$ should bewe will get something of the form $a+b \pi + c \alpha + d \beta$$a+b \pi$, where $\alpha$ and $\beta$ are the periods of $E$as desired. The classical way
Actually, it looks to state this is that everything can be expressed in terms ofme like we should just get complete elliptic integrals of the first, second and third kind$b \pi$. Maybe the integration by parts doesn't go as well as I hoped?
But I don't know why you aren't seeing those period termsObviously, andsomeone should actually work this out explicitly, but I don't know how to make this explicit. So I'm going to stop herethink it will be me.