Consider the following oscillatory integral $$ I(n):=\int_{-\pi}^\pi\int_{-\pi}^\pi e^{i n(x+y)}\frac {(1 - \cos(2x)) (1 - \cos(2y))} {2k - (\cos x + \cos y)}\ \mathrm{d}x\,\mathrm{d}y. $$ where $n\in\mathbb{N}$, $n>2$, $k\in\mathbb{R}$, $k>1$, and $i$ is the imaginary unit.
My question. Is it possible to characterize the asymptotic decay rate of $I(n)$ as a function of the integer parameter $n$?
The asymptotics for large $n$ is $$ J(n,\kappa) := \Big(\frac{2}{\pi}\Big)^2 \int_{-\pi}^\pi \int_{-\pi}^\pi \exp{(i\,n(x+y))}\frac{\sin^2x\,\sin^2y} {2\kappa - (\cos{x}+\cos{y}) }\, dx \,dy \sim$$ $$ \sim \frac{8}{\sqrt{\pi \kappa n}}(\kappa^2-1)^{7/4} (\kappa - \sqrt{\kappa^2-1})^{2n}\quad, \quad (\kappa>1)$$ The proof consists of 5 parts. The first part is to write the 2-dimensional fourier transform in a 1-dimensional form, albeit with special functions as part of the integrand. In Parts 2-5 asymptotic analysis will be undertaken. The proof is already long, and I'm not going to justify every switch of $\sum$ and $\int$, nor do some of the details that a rigorous saddle point analysis requires, like proving the tails decay sufficiently fast. I have checked the result numerically, and for $\kappa$ = 2.2, I get the following errors between the true value and the asymptotic approximation: n=50, 2.46%; n=100, 1.23%, n=200, 0.61%
Part 1: Show $$ J(n,\kappa) = \int_0^\infty \exp{(-2\kappa t)}\Big(2I_n(t) - I_{n-2}(t) - I_{n+2}(t) \Big)^2 dt $$ where $I_n(t)$ is the modified Bessel function.
With $r=1/(2\kappa)$,expand the integrand in a power series in $r$ to get $$ J(n,\kappa) =r \sum_{m=0}^\infty r^m \Big(\frac{2}{\pi}\Big)^2 \int_{-\pi}^\pi \int_{-\pi}^\pi \big(\cos{x}+\cos{y})^m \sin^2{x} \sin^2{y} \,e^{i\,n(x+y)} dx \, dy =$$ $$ = r \sum_{m=0}^\infty r^m m! \sum_{k=0}^m b_{m-k} \, b_k \, \text{ where }b_k=\frac{2}{k!\,\pi}\int_{-\pi}^\pi \cos^k y \sin^2 y \, e^{i\,n y} dy$$ By Cauchy product of power series, where $B(u)=\sum_k u^k\,b_k,$ $$J(n,\kappa) =r \sum_{m=0}^\infty r^m m! [u^m] B(u)^2 $$ where the 'coefficient of' operator has been used. By the Borel transform we get $$J(n,\kappa) = \int_0^\infty dt \exp{(-2\kappa t)} \sum_{m=0}^\infty t^m [u^m] B(u)^2 = \int_0^\infty \exp{(-2\kappa t)}B(t)^2 dt$$ We work on $B(t)$ now. $$ B(t) = \frac{2}{\pi} \int_{-\pi}^\pi dy \, e^{i\,n\,y} \sin^2 y \sum_{k=0}^\infty\frac{(t\, \cos{y})^k}{k!} $$ The interior summation is $\exp{(t \cos{y} )}$ which also has an expansion involving modified Bessel functions $$ \exp{(t \cos{y} )} = I_0(t)+2\sum_{j=1}^\infty I_j(t) \cos{j\,y} $$ We thus must consider the integral $$\frac{2}{\pi} \int_{-\pi}^\pi dy \, e^{i n \, y} \sin^2{y} \cos{j\,y} $$ Mathematica knows that $$\frac{2}{\pi}\int_{-\pi}^\pi e^{i n \, y} \sin^2{y} \, dy = \frac{-2}{\pi} \frac{4 \sin{(\pi z)}}{z(z+2)(z-2)} $$ which, for integer $z,$ is non-zero only if $z=0$ or $\pm 2.$ The penultimate formula says we have $z=n \pm j,$ and the formula under the 'Part 1' heading follows by keeping only the non-zero terms in the infinite sum over Bessel functions.
Part 2: a useful simplification
$$2I_n(u)-I_{n-2}(u)-I_{n+2}(u)=\frac{2}{u}\Big( (n+1)I_{n+1}(u)-(n-1)I_{n-1}(u) \Big) $$ This follows from applying the recursion $I_{n-1}(u)-I_{n+1}(u)=\frac{2n}{u}I_n(u)$ twice. The equation of Part 1 will thus consist of two BesselI's squared, and a cross term.
Part 3: Show $$\int_0^\infty e^{-2 \kappa\,u}\,I_n(u)^2 \frac{du}{u^2}= \frac{1}{n\,\pi(4n^2-1)} \int_0^1 \frac{x^{-n}}{\sqrt{x(1-x)}}\big(\kappa-\sqrt{\kappa^2-x}\big)^{2n} \big(\kappa + 2n\sqrt{\kappa^2 - x} \big)\,dx$$
Use Gradshteyn 6.592.6 and 6.611.4, respectively: $$\frac{1}{\pi} \int_0^1 \frac{dx}{\sqrt{x(1-x)}}I_{2n}(2u\sqrt{x}) = I_n(u)^2 \, , $$ $$ \int_0^\infty e^{-a\,x} I_{\nu}(b\,x)\,dx = b^{-\nu}\frac{\big(a-\sqrt{a^2-b^2}\big)^{\nu}}{\sqrt{a^2-b^2}} .$$
Then integrate twice with respect to the argument of the exponential ($a,$ to begin with).
Part 4: Asymptotics of 'cross term'
$$M_n:=\int_0^\infty e^{-2 \kappa\,u}\,I_{n-1}(u)\,I_{n+1}(u) \frac{du}{u^2}=
\int_0^\infty e^{-2 \kappa\,u}\,I_{n}(u)^2 \frac{du}{u^2} \big(1+\cal{O}(1/n)\big) $$
Use the series form for product of Bessel functions, Gradsteyn 8.442.1. Then
$$M_n=\int_0^\infty \frac{du}{u^2}e^{-2 \kappa\,u} \sum_{m=0}^\infty \frac{(u/2)^{2n+2m}}{m!} \frac{(2n+m+1)_m}{(n+m+1)!(n+m-1)!} =$$
$$=\int_0^\infty \frac{du}{u^2}e^{-2 \kappa\,u} \sum_{m=0}^\infty \frac{(u/2)^{2n+2m}}{m!} \frac{(2n+m+1)_m}{(n+m)!(n+m)!}\Big(1 - \frac{1}{n+m+1} \Big)$$
$$=\int_0^\infty \frac{du}{u^2}e^{-2 \kappa\,u} I_n(u)^2\big(1+\cal{O}(1/n)\big) $$
Part 5: Put pieces together
$$ J(n,\kappa) \sim 4\int_0^\infty \frac{du}{u^2}e^{-2 \kappa u}\Big( \big((n+1)I_{n+1}(u)\big)^2 + \big((n-1)I_{n-1}(u)\big)^2 -2(n+1)(n-1)I_{n}(u)^2\Big)$$ and using the result of Part 3,
$$\frac{\pi}{4}J(n,\kappa)\sim \int_0^1 \frac{dx\,x^{-n}}{\sqrt{x(1-x)}}\big( \kappa - \sqrt{\kappa^2-x}\big)^{2n} \cdot$$
$$\cdot \Big\{ \frac{n+1}{4(n+1)^2-1} \frac{\big( \kappa - \sqrt{\kappa^2-x}\big)^2}{x}\big(\kappa+2(n+1)\sqrt{\kappa^2-x}\big) - $$
$$ \frac{n-1}{4(n-1)^2-1} \frac{x}{\big( \kappa - \sqrt{\kappa^2-x}\big)^2}\big(\kappa+2(n-1)\sqrt{\kappa^2-x}\big) + $$
$$ \frac{2(n+1)(n-1)}{n(4n^2-1)} \big(\kappa+2n\sqrt{\kappa^2-x}\big) \Big\} $$
Expand curly brackets as $n \to \infty.$
$$ \Big\{...\Big\}=\frac{2}{x}\big(\kappa^2-x\big)^{3/2}+\kappa(\kappa^2/x-1)\frac{1}{n}+...$$ Since earlier we approximated the cross-term to order 1/n, we'll keep only the first term in this expansion. By letting $x \to u^2$ then
$$ \frac{\pi}{16} J(n,\kappa) \sim \int_0^1 \frac{du \, u}{\sqrt{(1-u)(1+u)}}\big(\kappa/{u}-\sqrt{ (\kappa/u)^2-1}\big)^{2n}\big((\kappa/u)^2-1)^{3/2} $$
Use the following transformation, good for positive functions $g,$
$$\int_0^1g(t)/\sqrt{1-t}=2\int_0^1 g(4u(1-u))\,du$$
This has the effect of putting the strong growth behavior at the endpoint of the domain as $u \to 1$ to the region near $u=1/2.$ In fact, the following integral is well-setup for a classic saddle point analysis:
$$ \frac{\pi}{128}J(n,\kappa) \sim \int_0^1 \frac{du \, u(1-u)}{\sqrt{1+4u(1-u)}} \Big(\big(\frac{\kappa}{4u(1-u)}\big)^2-1\Big)^{3/2}\cdot $$
$$\Big(\frac{\kappa}{4u(1-u)}-\sqrt{ \big(\frac{\kappa}{4u(1-u)}\big)^2-1}\Big)^{2n}$$
Define
$$h(u)=-\log\Big(\frac{\kappa}{4u(1-u)}-\sqrt{ \big(\frac{\kappa}{4u(1-u)}\big)^2-1}\,\Big) $$
We find $h'(u)=0 \implies u=1/2.$ A Taylor expansion about u=1/2 leads to
$$ h(u)=-\log\Big(\kappa-\sqrt{\kappa^2-1}\Big) +\frac{4\kappa}{\sqrt{\kappa^2-1}}(u-1/2)^2+\cal{O}((u-1/2)^4) $$
Define $g(u)$ as below and keep its first term in a Taylor expansion about u=1/2
$$ g(u) = \frac{ u(1-u)}{\sqrt{1+4u(1-u)}} \Big(\big(\frac{\kappa}{4u(1-u)}\big)^2-1\Big)^{3/2}=\frac{\big(\kappa^2-1)^{3/2}}{4\sqrt{2}}+\cal{O}((u-1/2)^2)$$
Then
$$ \frac{\pi}{128}J(n,\kappa) \sim \Big(\kappa-\sqrt{\kappa^2-1}\Big)^{2n}
\frac{\big(\kappa^2-1)^{3/2}}{4\sqrt{2}} \int_{-\infty}^\infty
\exp{\Big(\frac{-8\kappa \,n}{\sqrt{\kappa^2-1}}(u-1/2)^2 \Big)} du $$
Use of the Gaussian integral and algebra completes the proof.
The leading order exponential can be (almost) found by Paley-Wiener (which also applies to Fourier series).
Your integrand extends to a complex analytic function on the strip $$|\Im x|, |\Im y| < a = \cosh^{-1}(k) = \ln(k + \sqrt{k^2 - 1})$$ So the Payley-Wiener theorem says that (after performing the Fourier transform separately in $x$ and $y$) for such functions the Fourier series has the decay $$ I(n) \lesssim_{\epsilon} e^{-2an + n\epsilon} $$ for any $\epsilon > 0$. (The factor of 2 coming from the fact that you are looking at the Fourier series in $x$ times the Fourier series in $y$.) Plugging in you get
$$ I(n) \lesssim_\epsilon \left(\frac{1}{k + \sqrt{k^2 - 1}} \right)^{2n} \cdot e^{n\epsilon} $$
(Note that $(k + \sqrt{k^2 - 1})^{-1} = k - \sqrt{k^2 - 1}$.)
For the precise leading order you will probably need to do some contour integrals by hand.
