Uma integral difícil

Seja $\Phi$ a solução fundamental da equação do calor, dada por

$$\Phi(x,t)=\left\{\begin{align*}
\frac{1}{(4\pi t)^{n/2}}\mathrm{e}^{\frac{-|x|^2}{4t}},&\qquad x\in\mathbb{R}^n,\; t>0\\
0,& \qquad x\in\mathbb{R}^n,\; t<0
\end{align*}\right.$$

Fixados  $x\in\mathbb{R}^n$ e $t\in\mathbb{R}$, defina

$$E(x,t;r)=\left\{(y,s)\in\mathbb{R}^{n+1}\;\Big |\; s\leq 0,\;\Phi(x-y,t-s)\geq \frac{1}{r^n}\right\}.$$

Problema: Calcular a integral múltipla abaixo. $$\iint_{E(0,0;1)}\frac{|y|^2}{s^2}\;dyds.\tag{$*$}$$

Observação: Para o caso $n=1$, o problema se reduz a calcular a integral dupla $\iint_{\mathcal{R}}\frac{x^2}{y^2}\;dxdy$, onde $\mathcal{R}$ representa a região

definida por $\mathcal{R}=\left\{(x,y)\in\mathbb{R}^{2}\mid -\tfrac{1}{4\pi}\leq y\leq 0,\;x^2\leq 2y\ln(-4\pi y)\right\}.$

Origem do Problema: O livro Partial Differential Equations de L. C. Evans, na Seção 2.3.2 (que trata sobre a equação do calor) utiliza (no meio da demonstração do Teorema 3) o fato de que a integral $(*)$ vale exatamente $4$. Apesar deste resultado não ser óbvio e nem fácil de obter, o livro não traz sequer uma palavra sobre como efetuar o cômputo desta integral. Nesta postagem apresentaremos o referido cômputo, que resulta de uma discussão ocorrida em 2013 no fórum Mathematics Stack Exchange.

Solução do Problema:  Para $x\in\mathbb{R}^n$ e $t>0$, temos

$$\Phi(x,t)=\frac{1}{(4\pi t)^{n/2}}\mathrm{e}^{\frac{-|x|^2}{4t}}$$

Portanto, para $y\in\mathbb{R}^n$ e $s<0$,

$$\begin{align*}
\Phi(-y,-s)\geq 1\quad\Longleftrightarrow& \quad\frac{1}{(-4\pi s)^{n/2}}\mathrm{e}^{\frac{|y|^2}{4s}}\geq1\\\\
\Longleftrightarrow& \quad(-4\pi s)^{n/2}\leq\mathrm{e}^{\frac{|y|^2}{4s}}\\\\
\Longrightarrow& \quad(-4\pi s)^{n/2}\leq 1\\\\
\Longleftrightarrow& \quad s\geq -\frac{1}{4\pi}
\end{align*}$$
e
$$\begin{align*}
\Phi(-y,-s)\geq 1\quad\Longleftrightarrow& \quad\ln\left(\frac{1}{(-4\pi s)^{n/2}}\mathrm{e}^{\frac{|y|^2}{4s}}\right)\geq0\\\\
\Longleftrightarrow& \quad(-n/2)\ln\left(-4\pi s\right)+\frac{|y|^2}{4s}\geq0\\\\
\Longleftrightarrow& \quad|y|^2\leq 2ns\ln(-4\pi s)
\end{align*}$$
Isto implica que
$$\begin{align*}E(0,0;1)&=\{(y,s)\in\mathbb{R}^{n+1}\mid s\leq 0,\;\Phi(-y,-s)\geq 1\}\\\\
&=\left\{(y,s)\in\mathbb{R}^{n+1}\;\Big|\; -\frac{1}{4\pi}\leq s\leq 0,\;|y|^2\leq 2ns\ln(-4\pi s)\right\}\end{align*}$$
Assim,
$$\begin{align*}
\iint_{E(0,0;1)}\frac{|y|^2}{s^2}\;dyds&=\int_{-\frac{1}{4\pi}}^0\int_{|y|^2\leq2ns\ln(-4\pi s)}\frac{|y|^2}{s^2}\;dyds\\\\
&=\int_{-\frac{1}{4\pi}}^0\int_0^{\sqrt{2ns\ln(-4\pi s)}}\int_{\partial B(0,r)}\frac{|y|^2}{s^2}\;dS(y)drds \quad\text{(coordenadas polares)}\\\\
&=\int_{-\frac{1}{4\pi}}^0\int_0^{\sqrt{2ns\ln(-4\pi s)}}\int_{\partial B(0,1)}\frac{|rw|^2}{s^2}r^{n-1}\;dS(w)drds \quad\text{(mudança $w=y/r$)}\\\\
&=\int_{-\frac{1}{4\pi}}^0\int_0^{\sqrt{2ns\ln(-4\pi s)}}\int_{\partial B(0,1)}\frac{r^{n+1}}{s^2}\;dS(w)drds\\\\
&=\text{med}\big(\partial B(0,1)\big)\int_{-\frac{1}{4\pi}}^0\int_0^{\sqrt{2ns\ln(-4\pi s)}}\frac{r^{n+1}}{s^2}drds\\\\
&=\text{med}\big(\partial B(0,1)\big)\int_{-\frac{1}{4\pi}}^0\frac{(2ns\ln(-4\pi s))^{(n+2)/2}}{(n+2)s^2}ds\\\\
&=\frac{(2n)^{(n+2)/2}}{(n+2)}\text{med}\big(\partial B(0,1)\big)\int_{-\frac{1}{4\pi}}^0\frac{(s\ln(-4\pi s))^{(n+2)/2}}{s^2}ds
\end{align*}$$
Calculando esta última integral, obtemos
$$\begin{align*}
\int_{-\frac{1}{4\pi}}^0\frac{(s\ln(-4\pi s))^{(n+2)/2}}{s^2}ds&=\int_{1}^0\frac{\left((-\frac{x}{4\pi})\ln(x)\right)^{(n+2)/2}}{(-\frac{x}{4\pi})^2}\left(-\frac{1}{4\pi}\right)dx\quad\text{(mudança $x=-4\pi s$)}\\\\
&=\int_{0}^\infty\frac{\left(\frac{\mathrm{e}^{-z}}{4\pi}z\right)^{(n+2)/2}}{(-\frac{\mathrm{e}^{-z}}{4\pi})^2}\left(-\frac{1}{4\pi}\right)(-\mathrm{e}^{-z})\;dz\quad\text{(mudança $z=-\ln(x)$)}\\\\
&=\frac{1}{(4\pi)^{n/2}}\int_{0}^\infty\mathrm{e}^{-nz/2}z^{(n+2)/2}\;dz\\\\
&=\frac{1}{(4\pi)^{n/2}}\int_{0}^\infty\mathrm{e}^{-w}\left(\frac{2w}{n}\right)^{(n+2)/2}\frac{2}{n}\;dw\quad\text{(mudança $w=nz/2$)}\\\\
&=\frac{1}{(4\pi)^{n/2}}\left(\frac{2}{n}\right)^{n/2+2}\int_{0}^\infty\mathrm{e}^{-w}w^{n/2+1}\;dw\\\\
&=\frac{1}{2^{n}\pi^{n/2}2^{-n/2-2}n^{n/2+2}}\int_{0}^\infty\mathrm{e}^{-w}w^{n/2+2-1}\;dw\\\\
&=\frac{1}{\pi^{n/2}2^{n/2-2}n^{n/2+2}}\Gamma\left(2+\frac{n}{2}\right)\quad\text{(definição da função gama)}
\end{align*}$$
Deste modo,
$$\begin{align*}
\iint_{E(0,0;1)}\frac{|y|^2}{s^2}\;dyds&=\frac{(2n)^{(n+2)/2}}{(n+2)}\text{med}\big(\partial B(0,1)\big)\frac{1}{\pi^{n/2}2^{n/2-2}n^{n/2+2}}\Gamma\left(2+\frac{n}{2}\right)\\\\
&=\frac{8\;\text{med}\big(\partial B(0,1)\big)}{n(n+2)\pi^{n/2}}\Gamma\left(2+\frac{n}{2}\right)\\\\
&=\frac{8\;\text{med}\big(\partial B(0,1)\big)}{n(n+2)\pi^{n/2}}\left(1+\frac{n}{2}\right)\Gamma\left(1+\frac{n}{2}\right)\quad\text{(propriedade $\Gamma(t+1)=t\Gamma(t)$)}\\\\
&=\frac{4\;\text{med}\big(\partial B(0,1)\big)}{n\pi^{n/2}}\Gamma\left(1+\frac{n}{2}\right)\\\\
&=\frac{4\;\text{med}\big(\partial B(0,1)\big)}{n\pi^{n/2}}\frac{n}{2}\Gamma\left(\frac{n}{2}\right)\quad\text{(mesma propriedade anterior)}\\\\
&=4\ \text{med}\big(\partial B(0,1)\big)\frac{\Gamma\left(\tfrac{n}{2}\right)}{2\pi^{n/2}}\\\\
&=4\quad \text{(propriedade med$(\partial B(0,1))=\tfrac{2\pi^{n/2}}{\Gamma(n/2)}$)}
\end{align*}$$ 

Nota: A definição e as propriedades utilizadas da função gama podem ser encontradas nas páginas 7 e 8 da segunda edição do livro Introduction to partial differential equations de G. B. Folland.

Referências: aqui e aqui.