Integrate $\frac{R}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)} dx dy $

As in the title: let $m,n\in\mathbb{Z}$. Integrate:

$$ \frac{R}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)} dx dy $$ For me, it's quite a difficult problem. Any hints?

$R$ is a constant.

Here is an answer which partially identifies the closed form of the integral.

Also, since part 1 and 3 are too long, I tried to summarize them. You may skip details and jump to summaries.

1. Simplification

First, let $I(m,n)$ denote the integral. If $\mathcal{D}$ denotes the square with corners $(\pm 2\pi, \pm 2\pi)$, we have

\begin{align*} I(m, n) &= \frac{R}{8\pi^2} \iint_{\mathcal{D}} \frac{1-\cos(mx+ny)}{2-\cos x - \cos y} \, dxdy \\ &= \frac{R}{16\pi^2} \iint_{\mathcal{D}} \frac{1-\cos(mx+ny)}{1-\cos(\frac{x+y}{2})\cos(\frac{x-y}{2})} \, dxdy \\ &= \frac{R}{8\pi^2} \int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \frac{1-\cos((m+n)u+(m-n)v)}{1-\cos u \cos v} \, dudv \\ &= \frac{R}{2\pi^2} \int_{0}^{\pi}\int_{0}^{\pi} \frac{1-\cos((m+n)u)\cos((m-n)v)}{1-\cos u \cos v} \, dudv, \end{align*}

where $(u,v) = (\frac{x+y}{2}, \frac{x-y}{2})$. Now in view of the connection with the SRW Green function on $\Bbb{Z}^2$, it turns out that $R = 2$ is a natural choice. So we let

$$ A(m,n) = \frac{1}{\pi^2} \int_{0}^{\pi}\int_{0}^{\pi} \frac{1-\cos((m+n)u)\cos((m-n)v)}{1-\cos u \cos v} \, dudv \tag{$\diamond$} $$

so that $ I(m,n) = \frac{R}{2}A(m,n)$.

In order to determine the values of $A(m,n)$, note by symmetry that it suffices to consider the case $m \geq n \geq 0$. So we put $m = n+r$ with $n, r \geq 0$. Then using the formula

$$\int_{0}^{\pi} \frac{\cos(n\theta)}{1-r\cos\theta} \, d\theta = \frac{\pi}{\sqrt{1-r^2}}\left( \frac{1-\sqrt{1-r^2}}{r} \right)^n, \qquad |r| < 1 \wedge n \geq 0, $$

we can write

\begin{align*} A(n+r, n) &= \frac{1}{\pi} \int_{0}^{\pi} \frac{1}{\sin u} \left( 1 - \cos((2n+r)u) \left(\frac{1-\sin u}{\cos u}\right)^r \right) \, du \\ &= \frac{1}{\pi} \operatorname{Re} \int_{0}^{\pi} \frac{1}{\sin u} \left( 1 - e^{i(2n+r)u} \left(\frac{1-\sin u}{\cos u}\right)^r \right) \, du \\ &= \frac{2}{\pi} \operatorname{Re} \int_{C} \frac{1}{z^2 - 1} \left( 1 - z^{2n+r} \left(i \frac{z-i}{z+i}\right)^r \right) \, dz, \end{align*}

where $z = e^{i\theta}$ and $C$ is a contour from $1$ to $-1$ clockwise along the unit circle. Since the integrand has no singularity inside the unit disk, we may deform the contour to write

\begin{align*} A(n+r, n) &= \frac{2}{\pi} \operatorname{Re} \int_{-1}^{1} \frac{1}{1-x^2} \left( 1 - x^{2n+r} \left(i \frac{x-i}{x+i}\right)^r \right) \, dx \\ &= \color{red}{\frac{4}{\pi} \int_{0}^{1} \frac{1-x^{2n}}{1-x^2} \, dx + \frac{2}{\pi} \operatorname{Re} \int_{-1}^{1} x^{2n} \cdot \frac{1 - M_r(x)}{1-x^2} \, dx} \tag{1} \\ &= \color{red}{\frac{4}{\pi} \left( \int_{0}^{1} \frac{1-x^{2n}}{1-x^2} \, dx + \int_{0}^{1} x^{2n} \cdot \frac{1- \Re M_r(x)}{1-x^2} \, dx \right)}, \tag{2} \end{align*}

where $M_r(x)$ and $\Re M_r(x)$ are given by

$$ M_r(x) = x^{r} \left(i \frac{x-i}{x+i}\right)^r, \qquad \Re M_r(x) = \frac{M_r(x) + M_r(-x)}{2}. $$

When computing $\text{(1)}$, it turns out that the following quantities come in handy:

\begin{align*} a_n &= \int_{0}^{1} \frac{1-x^{2n}}{1-x^2} \, dx = \sum_{k=1}^{n} \frac{1}{2k-1}, \\ b_n &= \int_{0}^{1} \frac{x^{2n}}{1+x^2} \, dx = (-1)^n \left(\frac{\pi}{4} - \sum_{k=1}^{n} \frac{(-1)^{k-1}}{2k-1} \right). \end{align*}

Summary. If we define $A(m,n)$ by $(\diamond)$, then

• $I(m,n) = \frac{R}{2}A(m,n)$, and
• for each $n, r \geq 0$, $A(n+r,n)$ is given by the representation $\text{(1)}$: $$ A(m,n) = \frac{4}{\pi} \int_{0}^{1} \frac{1-x^{2n}}{1-x^2} \, dx + \frac{2}{\pi} \operatorname{Re} \int_{-1}^{1} x^{2n} \cdot \frac{1 - M_r(x)}{1-x^2} \, dx $$ where $M_r(x) = x^{r} \left(i \frac{x-i}{x+i}\right)^r$.

2. Computation for small $r$'s

In theory, the formula $\text{(1)}$ (or equivalently, $\text{(2)}$) allows us to compute $A(n+r,n)$ for each given $r$. Here are some examples:

$\boxed{\text{Case 1}} \ $ Consider the case $r = 0$. Then $\Re M_0(x) \equiv 1$ and hence by $\text{(2)}$,

$$ \color{blue}{A(n,n) = \frac{4}{\pi} a_n.} $$

$\boxed{\text{Case 2}} \ $ If $r = 1$, then we have

$$ \Re M_1(x) = \frac{2x^2}{x^2 +1} \qquad \Rightarrow \qquad \frac{1-\Re M_1(x)}{1-x^2} = \frac{1}{x^2 + 1}. $$

Consequently

$$ \color{blue}{A(n+1, n) = \frac{4}{\pi}(a_n + b_n).} $$

$\boxed{\text{Case 3}} \ $ If $r = 2$, then it is easy to check that

$$ \frac{1-\Re M_2(x)}{1-x^2} = \frac{1 + 4x^2 - x^4}{(x^2 + 1)^2} = -1 + \frac{4x^2}{(x^2+1)^2} + \frac{2}{x^2 + 1}. $$

Now the following simple recurrence relation is useful for our computation:

$$ \int_{0}^{1} \frac{x^{2n+2}}{(x^2+1)^{k+1}} \, dx = -\frac{1}{2k \cdot 2^k} + \frac{2n+1}{2k} \int_{0}^{1} \frac{x^{2n}}{(x^2+1)^k} \, dx. $$

From this it follows that

$$ \color{blue}{A(n+2,n) = \frac{4}{\pi}\left( a_n + (4n+4)b_n - \frac{2n+2}{2n+1} \right).} $$

3. General structure of $A(n+r,n)$

Although the aforementioned approach works for each given $r$, the issue is that this approach is still far from producing a good closed form that works for both $n$ and $r$ simultaneously. So in this part, we attempt to identify a rough structure of $A(n+r,n)$.

Notice that the singularities $x = \pm 1$ of $x^{2n} (1-M_r(x))/(1-x^2)$ are removable. Then this rational function has poles only at $-i$ and $\infty$, and hence we can apply the partial fraction decomposition to write as

$$ x^{2n} \frac{1-M_r(x)}{1-x^2} = \sum_{k=0}^{2n+r} c_{n,k} x^k + \sum_{k=1}^{r} \frac{d_{n,k}}{(x+i)^k} $$

for some $\{c_{n,k}\}_{k=0}^{2n+r-2}, \{d_{n,k}\}_{k=1}^{n} \subset \Bbb{Q}(i)$. This tells that

\begin{align*} A(n+r, n) &= \frac{4}{\pi} \bigg( a_n + \sum_{k=0}^{n+\lfloor r /2 \rfloor} \frac{1}{2n+2k+1} \operatorname{Re}(c_{n,k}) + \frac{\pi}{4} \operatorname{Im}(d_{n,1}) \\ &\hspace{4em} + \frac{1}{2}\operatorname{Re} \sum_{k=2}^{r} d_{n,k} \cdot \frac{(i-1)^{1-k} - (i+1)^{1-k}}{k-1} \bigg). \end{align*}

This proves that $A(n+r, n) \in \operatorname{Im}(d_{n,1}) + \pi \Bbb{Q}$. If we put $P_r(n) = (-1)^n \operatorname{Im}(d_{n,1})$, then this also proves that there exist rational numbers $P_r(n)$ and $Q_r(n)$ satisfying

$$ \color{red}{A(n+r,n) = \frac{4}{\pi}(a_n + P_r(n)b_n + R_r(n)).} $$

We remark that this decomposition is not arbitrary in the following sense

• $a_n$ is the leading order of $A(n+r,n)$ and the remaining terms vanishes as $n\to\infty$.
• $P_r(n)b_n$ produces the integer term, i.e., $A(n+r,n) \in (-1)^n P_r(n) + \frac{1}{\pi}\Bbb{Q}$.

Therefore the problem boils down to identifying $P_r(n)$ and $R_r(n)$.

At this moment, I was able to identify a simple formula for $P_r(n)$ and managed to prove it. Indeed, notice that $(x+i)^r x^{2n} (1 - M_r(x))/(1-x^2)$ is holomorphic near $x = -i$ and that $d_{n,1}$ is the coefficient of $(x+i)^{r-1}$ in the power series of this function near $x = -i$. So we may write

\begin{align*} P_r(n) &= \frac{(-1)^n}{(r-1)!} \operatorname{Im} \left. \frac{d^{r-1}}{dx^{r-1}} \right|_{x=-i} (x+i)^r x^{2n} \cdot \frac{1 - M_r(x)}{1-x^2} \\ &= \frac{(-1)^n}{(r-1)!} \operatorname{Im} \left. \frac{d^{r-1}}{dx^{r-1}} \right|_{x=-i} \frac{x^{2n+r} (1+ix)^r}{x^2-1}. \end{align*}

Using the general Leibniz rule, we easily compute that

$$ P_r(n) = \operatorname{Im} \sum_{k+l+m = r-1} \binom{2n+r}{k}\binom{r}{l} 2^k (-1+i)^{m+1}. $$

This is already good compared to the nebulous notation $d_{n,1}$, but we can go further. Notice that we can interpret this sum as the constant term of the following formal Laurent series

$$ \operatorname{Im} \left( (1+2x)^{2n+r} (1+x)^r \cdot \frac{-1+i}{1 - (-1+i)x} \cdot \frac{1}{x^{r-1}} \right) = \frac{(1+2x)^{2n+r}(1+x)^r}{(1+2x+2x^2)x^{r-1}}. \tag{*} $$

If we let $q = (1+2x)(1+x^{-1})$, we easily check that $\text{(*)}$ is written as

\begin{align*} \text{(*)} = (1+2x)^{2n} \cdot \frac{q^r}{q-1} &= \sum_{k=0}^{r-1}(1+2x)^{2n}q^k + \frac{(1+2x)^{2n}}{q-1} \\ &= \sum_{k=0}^{r-1}(1+2x)^{2n+k}(1+x^{-1})^k + \frac{x(1+2x)^{2n-1}}{1+x}. \end{align*}

Since the remainder term does not contribute to the constant term, it follows that

$$ \color{blue}{ P_r(n) = \sum_{k=0}^{r-1}\sum_{j=0}^{k} \binom{k}{j}\binom{k+2n}{j}2^j = \sum_{k=0}^{r-1} D(k,k+2n), } \tag{3} $$

where $D$ is the Delannoy numbers.

Summary. There exists rational numbers $P_r(n)$ and $R_r(n)$ such that $$ A(n+r,n) = \frac{4}{\pi}(a_n + P_r(n)b_n + R_r(n)). $$ This representation is meaningful in the following sense:

• $A(n+r, n) \sim \frac{4}{\pi} a_n$ in the sense that $P_r(n)b_n + R_r(n) \to 0$ as $n\to\infty$ for each $r \geq 0$.
• $P_r(n)$ is given by the formula $\text{(3)}$, and $(-1)^n P_r(n)$ is the rational part of the closed form of $A(n+r,n)$. I.e., $A(n+r,n) \in (-1)^n P_r(n) + \pi \Bbb{Q}$.

4. Conjecture

At this point I have no idea how to identify $R_r(n)$, but computation seems to suggest the following form for $R_r(n)$:

Conjecture. For each $r \geq 0$ there exists a polynomial $Q_r$ of non-negative coefficients such that for each $n \geq 0$, we have $Q_r(n) \in \Bbb{Z}$ and

$$ A(n+r,n) = \frac{4}{\pi} \left( a_n + P_r(n) b_n - Q_r(n)\prod_{j=1}^{\lfloor r/2\rfloor} \frac{1}{(2j-1)(2n+2j-1)} \right). \tag{4} $$

The previous section shows that $Q_0, Q_1, Q_2$ can be chosen as

$$Q_0(n) = 0, \qquad Q_1(n) = 0, \qquad Q_2(n) = 2n+2. $$

Computation using software produces the following results:

\begin{align*} Q_3(n) &= 4n^2 + 14n + 12 \\ Q_4(n) &= 32n^4 + 256n^3 + 776n^2+1060n+552 \\ Q_5(n) &= 32n^5 + 384n^4 + 1848n^3 + 4488n^2+5548n+2820 \\ & \vdots \end{align*}

(It is not true in general that $Q_r$ has integer coefficients, but still it seems that $Q_r$ maps $\Bbb{Z}$ into $\Bbb{Z}$.) Existence of a polynomial $Q_r(n)$ satisfying $\text{(4)}$ is not hard to prove using an appropriate partial fraction decomposition, but it seems not useful to establish an exact for $Q_r(n)$ or the conjecture.

Complete Closed form Solution

Defining:

$$A(m,n) \equiv \frac{1}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)} dx dy$$

and assuming that because of symmetry we may assume that $m \ge n \ge 0 $ I claim that the solution is:

$$\begin{align}A(m,n) &= \\ & (-1)^{n}\frac{1}{2} p_{(m-n-1)}\left(\frac{m+n}{2}\right)\tag{1a} \\ &+ \sum_{k=0}^{\left\lfloor \frac{(m+n-2)}{2} \right\rfloor}\frac{2\cos(\frac{\pi}{2}((m+n-2)- 2(k+\frac{m - n}{2})))}{\pi(2k+1)} p_{(m+n-2-2k)}(\frac{m - n}{2})\tag{1b} \\ &+ \sum_{k=0}^{\left\lfloor \frac{(m-n-2)}{2} \right\rfloor}\frac{2\cos(\frac{\pi}{2}((m-n-2)- 2(k+\frac{m + n}{2})))}{\pi(2k+1)} p_{(m-n-2-2k)}(\frac{m + n}{2})\tag{1c} \end{align}$$

In which $p_N(x)$ is an $N^{\text{th}}$ order polynomial, which can be written in recursive form as a weighted sum of the integrals of lower order polynomials:

$$p_N(x) = \cos(N\frac{\pi}{2}) + 4\sum_{j=1}^{\left\lceil \frac{N}{2} \right\rceil}\frac{\int p_{(N+1-2j)}(x) dx}{2j-1} $$

The 'closed form' of this polynomial is given by:

$$p_N(x) = \sum _{j=0}^{\left\lfloor \frac{N}{2} \right\rfloor} \frac{(4(x))^{(N-2j)}}{(N-2j)!} (-1)^j\sum_{k_1=0}^{j}\sum_{k_2=0}^{(j-k_1)} \cdots \sum_{k_{(N-2j)}=0}^{(j-\sum_{k\in\{k_1,\ldots k_{(N-2j)-1}\}}k)} \prod_{k\in\{k_1,\ldots k_{N-2j}\}}\frac{(-1)^k}{(2k+1)} $$

Argumentation

As argued in my other answer, the Laplacian of $A(m,n)$ given by:

$$ \begin{align} \mathbf{D}^2_{mn} \otimes A(m,n) &=\begin{bmatrix}0 & \frac 1 2 & 0\\ \frac 1 2 & -2 & \frac 1 2\\0 & \frac 1 2 & 0\end{bmatrix} \otimes A(m,n) \\ & = \delta(0,0) \end{align}$$

is zero everywhere, except at $(m,n) = (0,0)$ where it is 1.

Together with the following, trivial subset of solutions:

$$ A(m,n) = {\begin{cases} 0 & \text {if $(m,n)=(0,0)$} \\ 0.5 & \text {if $ m^2+n^2 = 1$} \\ \sum_{k=1}^m \frac {2} {\pi(2 k - 1)} & \text {if $ m=n$} \\ \end{cases}} $$

as already derived in Sangchul Lee's answer, this completely defines the whole set of solutions $A(m,n)$, which is argued below. Therefore, if we have a general solution that (1) yields the subset of trivial solutions, and (2) has a discrete Laplacian which is $\delta(0,0)$, then we have our solution.

Proof that (1a,b,c) yields the subset of trivial solutions.

$$\begin{align}A(0,0) &= (-1)^{0}\frac{1}{2} p_{(-1)}\left(\frac{0}{2}\right) + \sum_{k=0}^{-1}\text{etcetera} + \sum_{k=0}^{-1}\text{etcetera}\\ & = 0 \end{align}$$

$$\begin{align}A(1,0) &= \frac{1}{2} p_{0}\left(\frac{1}{2}\right) + \sum_{k=0}^{-1}\text{etcetera} + \sum_{k=0}^{-1}\text{etcetera}\\ & = \frac 1 2 \cos(0)\\ & = \frac 1 2 \end{align}$$

$$\begin{align}A(m,m) &= (-1)^{m}\frac{1}{2} p_{(-1)}(m) \\ &+ \sum_{k=0}^{m-1}\frac{2\cos(\frac{\pi}{2}((2m-2)- 2(k)))}{\pi(2k+1)} p_{(2m-2-2k)}(0) \\ &+ \sum_{k=0}^{-1}\text{etcetera} \\ \\&= 0 + \sum_{k=0}^{m-1}\frac{2\cos((m-1-k)\pi)}{\pi(2k+1)} p_{2(m-1-k)}(0) + 0 \\ &=\sum_{k=0}^{m-1}\frac{2\cos((m-1-k)\pi)}{\pi(2k+1)} \cos(2(m-1-k)\frac{\pi}{2}) \\ \\ &=\frac{2}{\pi}\sum_{k=0}^{m-1}\frac{1}{(2k+1)} \end{align}$$

$\blacksquare$

Proof that the Laplacian of (1a,b,c) is $\delta(0,0)$

From the trivial subset, it is easy to see that the Laplacian at $(m,n)=(0,0)$ equals $-2\times 0 + 4(\frac{1}{2}\times \frac{1}{2}) = 1$. Now for $m \ge n \ge 0$ and $m \neq 0$, we must show that:

$$4A(m,n) \overset{?}{=} A(m,n+1) + A(m,n-1) + A(m+1,n) + A(m-1,n)$$

Note that I claim that this property holds for all individual terms $\{(1a), (1b), (1c)\}$ of the given solution.

For now, I will only show this for the first term. So I must show that:

$$\begin{align} 4(-1)^{n}\frac{1}{2} p_{(m-n-1)}\left(\frac{m+n}{2}\right) \\ &\overset{?}{=} (-1)^{n+1}\frac{1}{2} p_{(m-n-2)}\left(\frac{m+n+1}{2}\right)\\&+ (-1)^{n-1}\frac{1}{2} p_{(m-n)}\left(\frac{m+n-1}{2}\right)\\&+ (-1)^{n}\frac{1}{2} p_{(m-n)}\left(\frac{m+1 +n}{2}\right)\\&+ (-1)^{n}\frac{1}{2} p_{(m-n-2)}\left(\frac{m-1+n}{2}\right) \\ \end{align}$$ $$\Leftrightarrow $$ $$\begin{align} 4 p_{(m-n-1)}\left(\frac{m+n}{2}\right) \\ &\overset{?}{=} \left(p_{(m-n)}\left(\frac{m+n}{2}+\frac{1}{2}\right) - p_{(m-n)}\left(\frac{m+n}{2}-\frac{1}{2}\right)\right) \\&- \\&\left(p_{(m-n-2)}\left(\frac{m+n}{2}+\frac{1}{2}\right) - p_{(m-n-2)}\left(\frac{m+n}{2}-\frac{1}{2}\right)\right) \end{align}$$

Now using the recursive formulation of the polynomial, and the notation $P_N(x) = \int p_N(x)dx$, we see that the constants ($\cos(N\frac{\pi}{2})$) cancel, so we must show that:

$$\begin{align} 4 p_{(m-n-1)}\left(\frac{m+n}{2}\right) \\ &\overset{?}{=} 4\sum_{j=1}^{\left\lceil \frac{m-n}{2} \right\rceil}\frac{1}{2j-1}\left.P_{(m-n+1-2j)}(x) \right|_{\left(\frac{m+n}{2}-\frac{1}{2}\right)}^{\left(\frac{m+n}{2}+\frac{1}{2}\right)} \\&- \\&4\sum_{j=1}^{\left\lceil \frac{m-n-2}{2} \right\rceil}\frac{1}{2j-1}\left.P_{(m-n-2+1-2j)}(x) \right|_{\left(\frac{m+n}{2}-\frac{1}{2}\right)}^{\left(\frac{m+n}{2}+\frac{1}{2}\right)} \end{align}$$

Combining the two definite integrals and dividing by four we get:

$$\begin{align} p_{(m-n-1)}\left(\frac{m+n}{2}\right) \\ &\overset{?}{=} \left.P_{(m-n-1)}(x)\right|_{\left(\frac{m+n}{2}-\frac{1}{2}\right)}^{\left(\frac{m+n}{2}+\frac{1}{2}\right)} \\ \\ &+ \sum_{j=2}^{\left\lceil \frac{m-n}{2} \right\rceil}\left.\left(\frac{-2}{(2j-1)(2j-3)}\right)P_{(m-n+1-2j)}(x)\right|_{\left(\frac{m+n}{2}-\frac{1}{2}\right)}^{\left(\frac{m+n}{2}+\frac{1}{2}\right)}\end{align}$$

Which rearranged says that the definite integral of this polynomial from $x_1 - \frac 1 2 $ to $x_1 + \frac 1 2 $ is the value of that polynomial at $x_1$ , plus the integral of its derivatives (using again the recusrive definition of the polynomials):

$$\begin{align} \left.P_{(m-n-1)}(x)\right|_{\left(x_1-\frac{1}{2}\right)}^{\left(x_1+\frac{1}{2}\right)} & = p_{(m-n-1)}(x_1) \\&+ \sum_{j=2}^{\left\lceil \frac{m-n}{2} \right\rceil} \left.\left(\frac{2}{(2j-1)(2j-3)}\right) P_{(m-n+1-2j)}(x)\right|_{\left(x_1-\frac{1}{2}\right)}^{\left(x_1+\frac{1}{2}\right)}\end{align}$$

$\blacksquare$

For the other two terms (1b) and (1c), the reasoning is similar.

Rationale behind the solution: Laplacian Minesweeping with a Faraday Shield

As already noted, from the subset of trivial solutions, combined with the harmonicity of $A(m,n)$, all other values can be calculated, like a game of minesweeper:

Sweeping all solutions from the trivial ones, and the Laplacian filter

Now it is illustrative to decompose the subset of trivial solutions as follows:

It turns out that only the multiples of $\frac{1}{2}$ are responsible for the Laplacian being 1 at $(0,0)$, and zero elsewhere. The diagonal elements involving the factors of $\frac{2}{\pi}$ will with the Laplacian filter all cancel to $0$.

The decomposition into the values $\{\frac{2}{\pi},\frac{8}{3\pi}, \frac{23}{15\pi},... \}$ can be continued by decomposing them in four diagonal, infinite halflines (see image below). Now the only way the Laplacian of a diagonal half-line can be zero, is that it ends in a perpendicular diagonal of alternating values (this can be easily verified):

Alternating diagonal acting like a 'Faraday Shield'

Behind this 'Faraday'-diagonal, we need a second alternating diagonal, to compensate for the first one. The main difference, is that its amplitude must increase linearly to both sides:

Second Faraday diagonal with linear increasing amplitude

In turn, this second diagonal must be compensated by a third diagonal, whose amplitude increases quadratically, and so on...

The diagonal alternating patterns with increasing amplitude (magnitude-scaled)

Once this idea is clear, it is easy to see that the needed polynomials describing the amplitudes of the oscillations, must be each others anti-derivatives. More precisely: each higher order polynomial is a weighted sum of the anti-derivatives of lower order polynomials as given above.

With a similar argument, we can also show that the four sides of the central diamond of halves are continued by infinite 'Faraday-lines'. Also these lead to higher order oscillations, which can be described by the same set of polynomials.

Numerical validation

I validated the integral numerically for several values of $(m,n)$ on a grid of $2^{13} \times 2^{13}$ gridpoints.

Image of the function $\frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)}$ as evaluated in the numerical validation

Matlab code

%% Set 0 <= n <= m
m = 3; n = 7

    % ensure that:  0 <= n <= m
    n = round(abs(n));
    m = round(abs(m));
    if n > m 
        a = n;
        n = m;
        m = a;
    end 

%% 1) Numerical integration (validation): 
N = 2^9;
[x,y] = meshgrid(2*pi*(-N/2:N/2-1)/N); 
f = 1/(2*pi)^2*(1 - cos(m*x + n*y))./(2 - cos(x) - cos(y)); 
f(N/2+1,N/2+1) =  1/(2*pi)^2*(n^2 + m^2); 
% compensate for discrete grid
f(N/2+1,N/2+1) = f(N/2+1,N/2+1)/2; 
figure; imagesc(f); axis equal tight;
Anm_numerical = sum(f(:))/N^2*4*pi^2


%% 2) Proposed analytical solution: 
% Initialize polynomial coefficients: 
PMAX = m+n - 1;
C = zeros(PMAX,PMAX); 
for k = 1:PMAX
    C(k,1) = cos((k-1)*pi/2);
    for kBack = k-1:-2:1
        C(k,2:end) = C(k,2:end) + 1/(k - kBack)*C(kBack,1:end-1);
    end
end

% 2a) irrational part: 
xVals = (m + [-n, +n])/2; % u and v
PVals = (m + [+n, -n]) - 2; % P_u and P_v

% Initialize integration sum: 
S = 0;
for ValsNr = 1:2
    x = xVals(ValsNr); % u and v
    P = PVals(ValsNr); % P_u and P_v

    % Initialize partial sum: 
    s = 0; 
    for k = 0:floor(P/2)
        pNr = P-2*k; % polynomial index (~=degree)
        value_pNr = 2/pi/(2*k+1); % 
        sign_pNr  = cos((pNr - 2*x)*pi/2); % sign of oscillating polynomial
        for powIndex = 0:floor(pNr/2)
            pow = pNr - 2*powIndex; % power of coefficient
            coef = (4*x)^pow/factorial(pow)*C(pNr+1, pow+1);
            s = s + coef*value_pNr*sign_pNr;
        end
    end
    S = S + s;
end

% 2b) rational part: 
x = (m + n)/2; % v
P = (m - n)-1; % P_v

% Initialize partial sum: 
s = 0; 
pNr = P; % polynomial index (~=degree)
value_pNr = 1/2; % 
sign_pNr  = (-1)^n; % sign of oscillating polynomial
for powIndex = 0:floor(pNr/2)
    pow = pNr - 2*powIndex; % power of coefficient
    coef = (4*x)^pow/factorial(pow)*C(pNr+1, pow+1);
    s = s + coef*value_pNr*sign_pNr;
end
S = S + s;

Anm_analytic = S
ErrorBetweenNumericalAnalytic = (Anm_analytic - Anm_numerical)

Approach

I will show that when $R=1$ the resulting discrete scalar distribution:

$$ A(m,n) = \frac{1}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)} dx dy $$

is harmonic everywhere except when $(m,n)=(0,0)$. (In fact, $A (m,n)$ is the inverse operator of the Discrete Laplacian). This allows us to find all integrals with recursion, once the simple forms are calculated.

The discrete Laplacian

Let:

$$ g(m,n) = \frac{1-\cos(mx+ny)}{2-(\cos x + \cos y)} $$

Then the Discrete Laplacian of $g(m,n)$ in the $(m,n)$ domain is:

$$ \begin{align} \mathbf{D}^2_{mn} \otimes g(m,n) &=\begin{bmatrix}0 & \frac 1 2 & 0\\ \frac 1 2 & -2 & \frac 1 2\\0 & \frac 1 2 & 0\end{bmatrix} \otimes g(m,n) \\ & = \cos(mx+ny) \end{align}$$

This is not to hard to proof: Wolfram alpha

Harmonicity

Now define:

$$ \begin{align} I(m,n) &\equiv \mathbf{D}^2_{mn} \otimes A (m,n) \\&=\mathbf{D}^2_{mn} \otimes \frac{1}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} g(m,n)dx dy \\ & = \frac{1}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} \cos(mx+ny) dx dy \\ &= \mathcal{F}(1)\\ &= {\begin{cases} 1 & \text {if $(m,n)=(0,0)$} \\ 0 & \text{otherwise} \\ \end{cases}} \end{align} $$

which is the identity operator for convolution. This means that for each pair of integers $(m,n) \neq (0,0)$ the integral is the average of the four integrals of its direct neighbors $\{(m-1,n), (m+1,n), (m,n-1), (m,n+1)\}$.

Startvalues

It is easy to see that for $(m,n) = (0,0)$ the integral is zero. For the four direct neigbors $m^2 + n^2 = 1$, the integral must be $\frac 12$, because the Laplacian at $(0,0)$ is one, and because of symmetry.

Together with the integral derived for $m=m$ in the answer of Sangchul we have enough initial values to recursively find the other integrals:

$$ \frac{1}{4 \pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} g(m,n) dx dy = {\begin{cases} 0 & \text {if $(m,n)=(0,0)$} \\ 0.5 & \text {if $ m^2+n^2 = 1$} \\ \sum_{k=1}^m \frac {2} {\pi(2 k - 1)} & \text {if $ m=n$} \\ \end{cases}} $$

Recursive recipe:

If $n \neq m$, assume that $n>m$ (symmetry) then:

$$A(m,n) = 4*A(m,n-1) - (A(m-1,n-1) + A(m+1,n-1) + A(m,n-2)) $$

$\blacksquare$

Just an addendum to Sangchul Lee's answer.

If we exploit the parity of the cosine function we get that the wanted integral equals:

$$ I(m,n) = \frac{R}{\pi^2}\int_{0}^{\pi}\int_{0}^{\pi}\frac{1-\cos(nx)\cos(my)}{2-(\cos x+\cos y)}\,dy\,dy $$ and by writing $\frac{1}{2-(\cos x+\cos y)}$ as $\int_{0}^{+\infty}\exp\left[\left(\cos x+\cos y-2\right)t\right]\,dt$ we get: $$ I(m,n) = R\int_{0}^{\infty}e^{-2t}\left(I_0(t)^2-I_n(t)\,I_m(t)\right)\,dt $$ where $I_k$ is a modified Bessel function of the first kind. Since the inverse Laplace transform of $e^{-2t}$ is just $\delta(s-2)$, the problem boils down to computing the Laplace transform of $I_0(t)^2$ and $I_n(t)\,I_m(t)$. The Laplace transform of $I_0(t)^2$ is related with the complete elliptic integral of the first kind, and has a singularity at $s=2$, where it behaves like $-\frac{\log|2-s|}{2\pi}$.

In a explicit way, from: $$ I_0(z) = \sum_{n\geq 0}\frac{z^{2n}}{4^n n!^2} $$ it follows that, for any $\tau>2$: $$ I_0(z)^2 = \sum_{n\geq 0}\frac{\binom{2n}{n}}{4^n n!^2}\,z^{2n},\qquad \int_{0}^{+\infty}I_0(z)^2 e^{-\tau z}\,dz = \sum_{n\geq 0}\frac{\binom{2n}{n}^2}{4^n \tau^{2n+1}}=\frac{2}{\pi\tau}\,K\left(\frac{2}{\tau}\right) $$ At least in principle, we may perform the same computation for $I_m(t)\,I_n(t)$, then compute the Laplace transform of $I_0(t)^2-I_m(t)\,I_n(t)$ a $\tau=2$ through De l'Hopital theorem.