A limit in a Feynman "proof" about Fermat's Theorem.

As perhaps some of you already know, Richard P. Feynman, the famous physicist tried a non-orthodox (in his usual way, I suppose) proof of the Fermat's Last Theorem.

He tried a probabilistic "proof" that no rigurous mathematician in the world would have accepted, but that shows his enormous creativity and insights. It is very well explained here: http://www.lbatalha.com/blog/feynman-on-fermats-last-theorem

I have followed the derivation of Mr. Luis Batalha, but at some point he challenges the reader to prove that the value:

$$c_{n} = \int_{1}^{\infty}\int_{1}^{\infty}(u^{n} + v^{n})^{-1 + \frac{1}{n}} \, du\,dv $$

when $n \to \infty $ is approximately $1/n$ ($c_{n} \approx 1/n$).

Well, Mr. Batalha says that for big $n$ tends to $1/n$, and I think we can say for $n \to \infty$.

I'm afraid I have tried to solve the limit but I am clueless.

Thank you and congratulations to Mr. Batalha for such an interesting post.

EDIT: Typo in the definition of $c_{n}$, the low limit of the integral is 1, not 0 as it was before. I'm afraid the typo is also present in the link.

EDIT 2: When I asked this question a few days ago, I couldn't imagine so many rich and fruitful comments and answers. Thank you to all. I am going to let the question open for a few days more if someone wants to continue adding solutions. I think all the comments have been superb. I think I am goingo to choose the answer by @tired. Thank you all.

I can give a heuristic derivation of the results found numerically by @RaymondManzoni and @Arentino. I will try to make things more rigorous the next days

First of all we observe the invariance of the integral under the transformation $x \leftrightarrow y$ (i relabel $u\rightarrow x, v\rightarrow y $). Therefore we can write by symmetry

$$ I_n=\color{blue}{2}\int_1^{\infty}dx\int_x^{\infty}dy(x^n+y^n)^{-1+1/n} $$

now setting $(y/x)^n=q$ we obtain

$$ I_n=\frac{\color{blue}{2}}{n}\int_1^{\infty}dx\frac{1}{x^{n-2}}\int_1^{\infty}dq(1+q)^{-1+1/n}q^{-1+1/n} $$

the integral over $x$ is trivial and yields $$ I_n=\frac{\color{blue}{2}}{n}\frac{1}{n-3}\int_1^{\infty}dq(1+q)^{-1+1/n}q^{-1+1/n} $$

the trick now is to recognize that the $n$-dependence of the remaining integral is very weak if $n$ gets big (to be more precise it decreases in powers of $1/n$) so to leading order we may just ignore it $$ I_n\sim\frac{\color{blue}{2}}{n}\frac{1}{n-3}\int_1^{\infty}dq(1+q)^{-1}q^{-1}+\mathcal{O}\left(\frac{1}{n^3}\right) $$

and therefore

$$ I_n\sim\frac{\color{blue}{2}}{n}\frac{1}{n-3}\log(2)+\mathcal{O}\left(\frac{1}{n^3}\right)\sim\frac{\color{blue}{2}\log(2)}{n^2}+\mathcal{O}\left(\frac{1}{n^3}\right) $$

in agreement with numerical evaluations, and FEYNMAN SEEMS TO BE WRONG ;)

Remark

By pushing the expansions one step further, i can also confirm Raymonds calculations up to $1/n^3$ so i'm sure now that my method is correct

In terms of polar coordinates, $(u^n+v^n)^{-1+1/n} \leq r^{-n+1} 2^{(-n/2+1)(-1+1/n)} = r^{-n+1} 2^{n/2-3/2+1/n}$. This is attained along $\theta=\pi/4$. This makes geometric sense: the $n$ norms with $n \geq 2$ are smaller than the $2$ norm, and the place where they are smaller by the largest amount is along the lines $y=\pm x$. We're dealing with a negative power of an $n$ norm, so for fixed $r$ it is biggest along those lines.

This looks bad, because of that growing $2^{n/2}$. But you can rewrite this as $2^{-3/2+1/n} r (r/\sqrt{2})^{-n}$. On $[1,\infty) \times [1,\infty)$ you have $r \geq \sqrt{2}$. Thus the integral is bounded by a constant times $\int_1^\infty r^{2-n} dr$. This bound is finite as soon as $n>3$ and behaves like $\frac{1}{n}$ as $n \to \infty$.

I'm not sure whether the original integral is finite when $n=3$. (Also, I neither claim nor expect that these estimates are at all tight.)

The integral can be computed with Mathematica and the result (for $n>3$) is:

$$ \begin{align} c_n=&{1\over n-3} \left({4^{-1/n}\over\sqrt\pi} \Gamma\left({1\over2} - {1\over n}\right) \Gamma\left({1} + {1\over n}\right) -{}_2F_1\left({1\over n}, {n-1\over n}; {n+1\over n}; -1\right)\\ +{1\over n-2}\ {}_2F_1\left({n-2\over n}, {n-1\over n}; {2n-2\over n}; -1\right) \right),\\ \end{align} $$ where ${}_2F_1$ is the hypergeometric function.

The first term $1/(n-3)$ is nothing but $\int_1^\infty r^{2-n} dr$, as found by Ian in his answer. Unfortunately, however, the term in parentheses does not tend to a positive constant for $n\to\infty$, but is rather proportional to $1/n$, because:

$$ \lim\limits_{n\to\infty}n^2c_n =2\ln2. $$ It follows that Batalha's claim is false.