Analytic continuation of double integral $\int_0^\infty \int_0^\infty f(x,y)^{-s} dx dy$

Question

Consider the following integral $$F(s) = \int_0^\infty \int_0^\infty [xy(x+1) (y+1) (y+2) (x+y+1) (2 x+y+2 )(x+y+2))]^{-s} dxdy$$

It converges and thus defines an analytic function when $1/4<\Re(s)<1$, it can be meromorphically continued to $\mathbb{C}$.

Question: it is true that $F(s)$ is analytic at each $s=0,-1,-2\cdots$? Is it true that $F(-n) \in \mathbb{Q}$?

After some guessing, I think the following is true: $$F(0) \stackrel{?}{=} \frac{25}{32}$$ I am also interested in a formula that calculates $F(-n)$. A technique that generalizes to similar problems like $$\int_0^\infty \int_0^\infty [xy(x+a_1) (y+a_2) (y+a_3) (x+y+a_4) (2 x+y+a_5 )(x+y+a_6))]^{-s} dxdy\qquad a_i>0$$ would be very nice.

---

Remark: an approach that doesn't work well is trying to "explicitly evaluate" $F(s)$ in terms of known function, one doesn't expect it to have a simple expression.

Answer 1

Here is a proof that $F(0) = \frac{25}{32}$. The method below can also prove that $F(-n) \in \mathbb{Q}$ for all $n \in \mathbb{N}_0$ and, in principle, yield a systematic formula for computing its values. Using Mathematica, if I implemented everything correct, we get:

$$ \begin{array}{|c|cccc|} n & 0 & 1 & 2 & 3 \\ \hline F(-n) & \dfrac{25}{36} & \dfrac{5879}{1075200} & \dfrac{87085}{290594304} & \dfrac{44523917}{5215826739200} \end{array} $$

Since my answer is already long enough, let me leave the proof of $F(-n) \in \mathbb{Q}$ to the future me. (The idea is quite similar to that in the computation of $F(0)$, though. It just becomes more complicated.)

1. Decomposition

Assume for a moment that $\frac{1}{5}<\Re(s)<1$. We fix $y > 0$ and consider the function $$ x \mapsto \frac{1}{(-x)^s (-x-1)^s (2x + y + 2)^s (x+y+1)^s (x+y+2)^s}. $$

The branch cut of this function is $(-\infty, -1-\frac{1}{y}] \cup [-1, \infty)$. Also, for $s$ in the regime considered, none of the singularities, including $\infty$ harm the integrability of the function. Consequently, by limiting argument, its integral along the following contour vanishes:

(Consider a right-half keyhole contour and let the radius tend to infinity.) This particular choice of contour is based on the symmetry of zeros. Later, we will see that this symmetry greatly helps simplify the representation of the integral. Now, the integral along each colored contour can be simplified using the standard tricks, yielding

$$ \begin{align*} &\int_{0}^{\infty} \frac{\mathrm{d}x}{x^s (x+1)^s (2x + y + 2)^s (x+y+1)^s (x+y+2)^s} \\ &= \color{red}{\frac{1}{2i\sin(2\pi s)} \int_{\mathcal{H}^-} \frac{\mathrm{d}x}{(-x)^s (-x-1)^s (2x + y + 2)^s (x+y+1)^s (x+y+2)^s}} \\ &= \color{blue}{ - \frac{\sin(\pi s)}{\sin(2\pi s)} \int_{0}^{1} \frac{\mathrm{d}x}{x^s (1-x)^s (y + 2-2x)^s (y+1-x)^s (y+2-x)^s} } \\ &\qquad + \color{magenta}{ \frac{1}{2\sin(2\pi s)} \int_{-\infty}^{\infty} \frac{\mathrm{d}x}{(1+\frac{y}{2}-ix)^s (\frac{y}{2}-ix)^s (2ix)^s (\frac{y}{2}+ix)^s (1+\frac{y}{2}+ix)^s} } \\ &= - \frac{\sin(\pi s)}{\sin(2\pi s)} \int_{0}^{1} \frac{\mathrm{d}x}{x^s (1-x)^s (y + 2x)^s (y+x)^s (y+1+x)^s} \tag{1a}\label{e:1a}\\ &\qquad + \frac{1}{2\sin(2\pi s)} \int_{-\infty}^{\infty} \frac{\mathrm{d}x}{((1+\frac{y}{2})^2+x^2)^s ((\frac{y}{2})^2+x^2)^s (2ix)^s}. \tag{1b}\label{e:1b} \end{align*} $$ We investigate the integrals $\eqref{e:1a}$ and $\eqref{e:1b}$ separately. It turns out that $\eqref{e:1a}$ is easy to handle with, as it easily extends to a meromorphic function. As a challenge-lover, let us study $\eqref{e:1b}$ first.

2. Contribution from $\eqref{e:1b}$

The integral term in $\eqref{e:1b}$ further simplifies to $$ \begin{align*} &\int_{-\infty}^{\infty} \frac{\mathrm{d}x}{((1+\frac{y}{2})^2+x^2)^s ((\frac{y}{2})^2+x^2)^s (2ix)^s} \\ &= \frac{1}{\Gamma(s)^3} \int_{0}^{\infty} \mathrm{d}t_1 \int_{0}^{\infty} \mathrm{d}t_2 \int_{0}^{\infty} \mathrm{d}t_3 \, (t_1 t_2 t_3)^{s-1} \\ &\qquad \times \int_{-\infty}^{\infty} \mathrm{d}x \, \exp\left( -\left(\left(1+\frac{y}{2}\right)^2 + x^2\right)t_1 -\left(\left(\frac{y}{2}\right)^2 + x^2\right)t_2 - 2ixt_3 \right) \\ &= \frac{1}{\Gamma(s)^3} \int_{0}^{\infty} \mathrm{d}t_1 \int_{0}^{\infty} \mathrm{d}t_2 \int_{0}^{\infty} \mathrm{d}t_3 \, (t_1 t_2 t_3)^{s-1} \\ &\qquad \times \sqrt{\frac{\pi}{t_1 + t_2}} \exp\left( - \frac{t_3^2}{t_1 + t_2} -\left(1+\frac{y}{2}\right)^2t_1 -\left(\frac{y}{2}\right)^2t_2 \right) \\ &= \frac{\sqrt{\pi}\Gamma(\frac{s}{2})}{2\Gamma(s)^3} \int_{0}^{\infty} \mathrm{d}t_1 \int_{0}^{\infty} \mathrm{d}t_2 \, (t_1 t_2)^{s-1} (t_1 + t_2)^{(s-1)/2} \exp\left( -\left(1+\frac{y}{2}\right)^2t_1 -\left(\frac{y}{2}\right)^2t_2 \right) \\ &= \frac{\sqrt{\pi}\Gamma(\frac{s}{2})}{2\Gamma(s)^3} \int_{0}^{1} \mathrm{d}p \int_{0}^{\infty} \mathrm{d}r \, r^{(5s-3)/2} (pq)^{s-1} \exp\left( -\left( \left(1+\frac{y}{2}\right)^2 p +\left(\frac{y}{2}\right)^2 q\right) r \right) \\ &= \frac{\sqrt{\pi}\Gamma(\frac{s}{2})\Gamma(\frac{5s-1}{2})}{2\Gamma(s)^3} \int_{0}^{1} \mathrm{d}p \, (pq)^{s-1} \left( \left(1+\frac{y}{2}\right)^2 p +\left(\frac{y}{2}\right)^2 q\right)^{\frac{1-5s}{2}} \\ &= \frac{\sqrt{\pi}\Gamma(\frac{s}{2})\Gamma(\frac{5s-1}{2})}{2\Gamma(s)^3} \int_{0}^{1} \mathrm{d}p \, (pq)^{s-1} \left( (1+y)p + \frac{y^2}{4} \right)^{\frac{1-5s}{2}} \end{align*} $$ The last integral defines a meromorphic function at least for $0 < \Re(s) < \frac{1}{3}$ based on the worst-case $y=0$. Substituting $p=\cos^2(\theta/2)=\frac{r+1}{2}$, we obtain two different representations of $\eqref{e:1b}$: $$ \begin{align*} \eqref{e:1b} &= \frac{1}{2\sin(2\pi s)}\frac{2^{3s-1}\sqrt{\pi}\Gamma(\frac{s}{2})\Gamma(\frac{5s-1}{2})}{\Gamma(s)^3} \int_{0}^{\pi} \mathrm{d}\theta \, (\sin\theta)^{2s-1} (y + 1 - e^{i\theta})^{\frac{1-5s}{2}} (y + 1 - e^{-i\theta})^{\frac{1-5s}{2}} \\ &= \frac{2^{2s-1}\Gamma(2s)\Gamma(1-2s)\Gamma(\frac{5s-1}{2})}{\Gamma(\frac{s+1}{2})\Gamma(s)^2} \int_{0}^{\pi} \mathrm{d}\theta \, (\sin\theta)^{2s-1} (y + 1 - e^{i\theta})^{\frac{1-5s}{2}} (y + 1 - e^{-i\theta})^{\frac{1-5s}{2}} \\ &= \frac{2^{2s-1}\Gamma(2s)\Gamma(1-2s)\Gamma(\frac{5s-1}{2})}{\Gamma(\frac{s+1}{2})\Gamma(s)^2} \int_{-1}^{1} \mathrm{d}r \, (1 - r^2)^{s-1} ((y+1)^2 - 2(y+1)r + 1)^{\frac{1-5s}{2}} \end{align*} $$

Below is a Mathematica code for numerically testing this equality:

s = 2/3;
y = 1/2;
(* Integral (1b) *)
1/(2 Sin[2 Pi  s])
  NIntegrate[
  1/(((1 + y/2)^2 + x^2)^s ((y/2)^2 + x^2)^s (2 I  x)^
   s), {x, -Infinity, Infinity}, WorkingPrecision -> 20]
(* Simplified version *)
((2^(2 s - 1) Gamma[2 s] Gamma[1 - 2 s] Gamma[(5 s - 1)/2])/(
  Gamma[(s + 1)/2] Gamma[s]^2)) NIntegrate[Sin[t]^(
  2 s - 1)/((y + 1 - Exp[I t])^((5 s - 1)/2) (y + 1 - Exp[-I  t])^((
   5 s - 1)/2)), {t, 0, Pi}, WorkingPrecision -> 20]

Now we integrate both sides of the above equality and apply the substitution $y \mapsto y^{-1} - 1$. The resulting integral can be written as $$ \begin{align*} \int_{0}^{\infty} \mathrm{d}y \, \frac{\eqref{e:1b}}{y^s(y+1)^s(y+2)^s} &= \frac{2^{2s-2}\Gamma(2s+1)\Gamma(1-2s)\Gamma(\frac{5s-1}{2})}{\Gamma(s+1)\Gamma(\frac{s+1}{2})} B_1(s), \tag{2.1}\label{e:2.1} \end{align*} $$ where the function $B_\rho(s)$ is defined for $\rho \in [0, 1]$ by $$ B_\rho(s) = \frac{1}{\Gamma(s)} \int_{-1}^{1} \mathrm{d}r \, (1 - r^2)^{s-1} \int_{0}^{\rho} \frac{y^{8s-3} \, \mathrm{d}y}{(1-y^2)^s (1 - 2ry + y^2)^{\frac{5s-1}{2}}}. \tag{2.2}\label{e:2.2} $$ The inclusion of the factor $\frac{1}{\Gamma(s)}$ is crucial, as we expect that $B(s)$ extends to a meromorphic function which has no pole at $s = 0$. Note also that the defining expression for $B(s)$ converges at least for $\frac{1}{4} < \Re(s) < \frac{1}{3}$.

2.1 Meromorphic Extension

This section is long, so you can only read the following summary and skip to the next section.

Summary.

• For $\rho < 1$, $B_\rho(s)$ extends to a meromorphic function on $\mathbb{C}$. We use the same symbol for its extension.
• $B_{\rho}(s)$ converges locally uniformly as $\rho \to 1^-$ for $\Re(s) < \frac{1}{3}$.
• The limiting function coincides with $B_1(s)$ on an open set. We also denote this extension by $B_1(s)$.

Our goal is to extend this to a meromorphic function on all of $\mathbb{C}$. For this, we consider the figure-eight curve $\gamma : [0, 1] \to \mathbb{C}$ winding $-1$ and $1$ as follows:

Now fix $\rho \in (0, 1)$, and choose $\varepsilon > 0$ and the curve $\gamma$ so that the following conditions are satisfied:

1. the zeros of $1 - 2zy + y^2$ for $z$ in the $\varepsilon$-neighborhood of $[-1, 1]$ are confined in a $(1-\rho)$-neighborhood of the unit circle. This is possible, because if $|y| \leq \rho$, then $1 - 2zy + y^2 = 0$ implies $z = \frac{1}{2}(y + y^{-1})$. It can be shown that the set of all such values are at positive distance from $[-1, 1]$.
2. $\gamma$ satisfies $\gamma(0) = \gamma(1) = 0$, and
3. $\gamma$ winds $1$ CCW first and then winds $-1$ CW.

Then we define a version of $\log(z^2 - 1)$ as follows: $$ L(t) = -\pi i + \int_{\gamma|_{[0, t]}} \left( \frac{1}{z + 1} + \frac{1}{z - 1} \right) \, \mathrm{d}z. $$ In terms of the principal complex logarithm, the following relation holds: $$ L(t) = \begin{cases} \log(1-\gamma(t)^2) - \pi i, & \text{along } \searrow \\ \log(1-\gamma(t)^2) + \pi i, & \text{along } \swarrow \end{cases} $$ (This is due to the branch cut $(-\infty, -1]\cup[1, \infty)$. Hence, $L(t)$ can be thought of as a continuous version of $\log(1-z^2)$ along $\gamma$, smoothly climing up and down of the logarithm "parking garage".) Under this setting, we have:

Lemma 1. Let $f(z)$ be analytic near $[-1, 1]$. If $\gamma$ is chosen to be close enough to $[-1, 1]$, then $$ \int_{-1}^{1} (1-x^2)^{s-1} f(x) \, \mathrm{d}x = \frac{1}{2i \sin(\pi s)} \int_{\gamma(t)} e^{(s-1)L(t)} f(z) \, \mathrm{d}z, $$ where we abuse the notation so that $\int_{\gamma(t)} f(z, t) \, \mathrm{d}z$ stands for $\int_{[0,1]} f(\gamma(t), t) \gamma'(t) \, \mathrm{d}t$.

Proof. Note that continuously deforming $\gamma$ in $\mathbb{C}\setminus\{\pm 1\}$ does not alter the value of the right-hand side. So by letting $\gamma$ to the union of two line segments, one from $-1+0^+i$ to $1-0^+i$ and the other from $1+0^+i$ to $-1 - 0^+i$, it follows that $$ \begin{align*} &\int_{\gamma(t)} e^{(s-1)L(t)} f(z) \, \mathrm{d}z \\ &= \int_{-1}^{1} e^{(s-1)(\log(1-x^2) - \pi i)} f(x) \, \mathrm{d}x + \int_{1}^{-1} e^{(s-1)(\log(1-x^2) + \pi i)} f(x) \, \mathrm{d}x. \end{align*} $$ Simplifying the right hand side, the desired conclusion follows.

Using this lemma, we can rewrite $B_{\rho}(s)$ as $$ B_\rho(s) = \frac{1}{2i \sin(\pi s)\Gamma(s)} \int_{\gamma(t)} \mathrm{d}z \, e^{(s-1)L(t)} \int_{0}^{\rho} \frac{y^{8s-3} \, \mathrm{d}y}{(1-y^2)^s (1 - 2zy + y^2)^{\frac{5s-1}{2}}}. \tag{2.3}\label{e:2.3} $$ This representation does not depend on the curve $\gamma$ so long as it satisfies all the technical conditions listed above. We further manipulate the integral by invoking:

Lemma 2. Let $f(z)$ be analytic near $|z| \leq \rho$. If $C(\rho)$ is the circle $|y| = \rho$ traced from $\rho e^{-i\pi}$ to $\rho e^{i\pi}$ in the CCW direction, then $$ \int_{0}^{\rho} x^{s-1} f(x) \, \mathrm{d}x = \frac{1}{2i \sin(\pi s)} \int_{C(\rho)} z^{s-1} f(-z) \, \mathrm{d}z. $$ Proof. Using the similar trick as before, we deform $C(\rho)$ to the union of two line segments, $\rho e^{-i \pi} \to 0 \to \rho e^{i\pi}$. (This notational abuse of distinguishing $e^{-i \pi}$ and $e^{i\pi}$ is most elegantly formalized when we consider the logarithmic Riemann surface, the covering space associated with $\exp$. Less elegantly, think of $e^{\pm i\pi}$ as the limit of $e^{\pm i(\pi - 0^+)}$, respectively.) $$ \begin{align*} \int_{C(\rho)} z^{s-1} f(-z) \, \mathrm{d}z &= \int_{0}^{\rho e^{i\pi}} z^{s-1}f(-z) \, \mathrm{d}z - \int_{0}^{\rho e^{-i\pi}} z^{s-1}f(-z) \, \mathrm{d}z \\ &= e^{i\pi s} \int_{0}^{\rho} x^{s-1}f(x) \, \mathrm{d}x - e^{-i\pi s} \int_{0}^{\rho} x^{s-1}f(x) \, \mathrm{d}x, \end{align*} $$ from which the desired conclusion follows.

Applying this lemma, $B_\rho(s)$ as $$ \begin{align*} B_\rho(s) &=- \frac{\Gamma(1-s)}{4\pi \sin(8\pi s)} \int_{\gamma(t)} \mathrm{d}z \, e^{(s-1)L(t)} \int_{C(\rho)} \frac{y^{8s-3} \, \mathrm{d}y}{(1-y^2)^s (1 + 2zy + y^2)^{\frac{5s-1}{2}}}. \tag{2.4}\label{e:2.4} \end{align*} $$ Here,. Note that the integrand has slightly changed, because we utilized the substitutions $y \mapsto e^{\pm i\pi}y$ in order to obtain $\eqref{e:2.4}$. The upshot of this calculation is that $\eqref{e:2.4}$ defines a meromorphic function on all of $\mathbb{C}$. This is because of all of $y$, $1-y^2$, and $1 + 2zy + y^2$ are bounded away from $0$ uniformly in $t \in [0, 1]$ and $y \in C(\rho)$.

We remark that $\eqref{e:2.4}$ is obtained only for each fixed $\rho \in (0, 1)$. The depenende of the choice of $\gamma$ on $\rho$ cannot be lifted, because the locus of zeros of $1+2zy + y^2$ gets closer to the unit circle as $\rho$ tends to $1$, essentially squeezing out the region available for $\gamma$. We cannot let zeros cross $\gamma$, since it would induce extra integral terms along a branch of $(1 + 2zy + y^2)^{\frac{5s-1}{2}}$, invalidating the representation.

However, $\eqref{e:2.4}$ is still good enough for our purpose. To this end, we show that $B_{\rho}(s)$ converges locally uniformly as $\rho \to 1^-$ on $\Re(s) < \frac{1}{3}$. To this end, let $0 < \rho_1 < \rho_2 < 1$. Then by $\eqref{e:2.3}$, $$ \begin{align*} B_{\rho_2}(s) - B_{\rho_1}(s) &= \frac{1}{2i \sin(\pi s)\Gamma(s)} \int_{\gamma(t)} \mathrm{d}z \, e^{(s-1)L(t)} \int_{\rho_1}^{\rho_2} \frac{y^{8s-3} \, \mathrm{d}y}{(1-y^2)^s (1 - 2zy + y^2)^{\frac{5s-1}{2}}}. \end{align*} $$ Now write $$ 1 - 2zy + y^2 = (1 + y)^2 \left(1 - \frac{2y}{(1+y)^2}(1+z) \right) $$ and note that $\frac{2y}{(1+y)^2} \leq \frac{2\rho_2}{(\rho_2 + 1)^2} < \frac{1}{2}$. This then implies $\left|\frac{2y}{(1+y)^2}(1+z)\right|$ can be bounded away from $1$ uniformly in $z \in \gamma([0, 1])$ and $y \in [\rho_1, \rho_2]$. Writing $p = \frac{5s-1}{2}$ for convenience, this allows us to expand $B_{\rho_2}(s) - B_{\rho_1}(s)$ as: $$ \begin{align*} B_{\rho_2}(s) - B_{\rho_1}(s) &= \frac{1}{2i \sin(\pi s)\Gamma(s)} \int_{\gamma(t)} \mathrm{d}z \, e^{(s-1)L(t)} \int_{\rho_1}^{\rho_2} \frac{y^{8s-3} \, \mathrm{d}y}{(1-y)^s(1+y)^{s+2p}} \\ &\hspace{5em} \times \sum_{n=0}^{\infty} (-1)^n \binom{-p}{n} \frac{(1+z)^n}{2^n} \left(\frac{4y}{(1+y)^2}\right)^n \end{align*} $$ By noting that $\left|\binom{-p}{n}\right| \sim \frac{n^{p-1}}{|\Gamma(p)|}$, the entire expression converges absolutely, hence we can rearrange the order of integrals and sums as we wish. We also observe that, applying Lemma 1 back, $$ \begin{align*} \frac{1}{2i \sin(\pi s)\Gamma(s)} \int_{\gamma(t)} \mathrm{d}z \, e^{(s-1)L(t)} \frac{(1+z)^n}{2^n \Gamma(s)} &= \int_{-1}^{1} \mathrm{d}r \, (1 - r^2)^{s-1} \frac{(1+r)^n}{2^n} \\ &= 2^{2s-1} \frac{\Gamma(n+s)}{\Gamma(n+2s)}. \end{align*} $$ Although the intermediate integral converges only for $\Re(s) > 1$, both the initial and final quantities are meromorphic on $\mathbb{C}$ in $s$, and so, the equality always holds accordingly. Then $$ \begin{align*} B_{\rho_2}(s) - B_{\rho_1}(s) &= \sum_{n=0}^{\infty} (-1)^n \binom{-p}{n} \int_{\rho_1}^{\rho_2} \frac{y^{8s-3} \, \mathrm{d}y}{(1-y)^s(1+y)^{s+2p}} \left(\frac{4y}{(1+y)^2}\right)^n \\ &\hspace{5em} \times \left[ \frac{1}{2i \sin(\pi s)} \int_{\gamma(t)} \mathrm{d}z \, e^{(s-1)L(t)} \frac{(1+z)^n}{2^n} \right] \\ &= \frac{2^{2s-1}}{\Gamma(p)} \sum_{n=0}^{\infty} \frac{\Gamma(n+p)}{n!} \frac{\Gamma(n+s)}{\Gamma(n+2s)} \\ &\hspace{5em} \times \int_{\rho_1}^{\rho_2} \frac{y^{8s-3} \, \mathrm{d}y}{(1-y)^s(1+y)^{s+2p}} \left(\frac{4y}{(1+y)^2}\right)^n. \end{align*} $$ Taking absolute value, $$ \begin{align*} \left| B_{\rho_2}(s) - B_{\rho_1}(s) \right| &\leq \left| \frac{2^{2s-1}}{\Gamma(p)} \right| \biggl( \int_{\rho_1}^{\rho_2} \mathrm{d}y \, \left| \frac{y^{8s-3}}{(1-y)^s(1+y)^{s+2p}} \right| \biggr) \sum_{n=0}^{\infty} \left| \frac{\Gamma(n+p)}{n!} \frac{\Gamma(n+s)}{\Gamma(n+2s)} \right| \end{align*} $$ By noting that $\frac{\Gamma(n+p)}{n!} \frac{\Gamma(n+s)}{\Gamma(n+2s)} \sim n^{\frac{3}{2}(s-1)}$ locally uniformly in $s$, the above bound converges to $0$ as $\rho_1, \rho_2 \to 1^-$ uniformly in the region $\Re(s) < \frac{1}{3}$ minus the poles of $\Gamma(s)$. The limit function is a meromorphic function at least for $\Re(s) < \frac{1}{3}$. However, we already know that $B_{\rho}(s) \to B_1(s)$ as $\rho \to 0^+$ locally uniform in $\Re(s) \in (\frac{1}{4}, \frac{1}{3})$. Hence we can unambiguously write $B_1(s)$ for the limiting function that extends the definition $\eqref{e:2.2}$ for $\rho = 1$.

2.2. Computing $B_1(0)$

In order to compute the value of $B_1(0)$, we compute $B_\rho(0)$ using $\eqref{e:2.4}$ and take limit as $\rho \to 1^-$. To this end, let $G_{\rho}(s, z)$ denote the innermost integral in $\eqref{e:2.4}$: $$ \begin{align*} G_\rho(s, z) := \int_{C(\rho)} \frac{y^{8s-3} \, \mathrm{d}y}{(1-y^2)^s (1 + 2zy + y^2)^{\frac{5s-1}{2}}} \end{align*} $$ When $s = 0$, this integral can be explicitly computed using residue: $$ \begin{align*} G_\rho(0, z) := \int_{C(\rho)} \mathrm{d}y \, \frac{\sqrt{1 + 2zy + y^2}}{y^3} = i \pi (1 - z^2) = -i\pi e^{L(z)}. \end{align*} $$ Consequently, the representation $\eqref{e:2.4}$ with $G_\rho(s, z)$ replaced by $G_\rho(0, z)$ is computed as $$ \begin{align*} &- \frac{\Gamma(1-s)}{4\pi \sin(8\pi s)} \int_{\gamma(t)} \mathrm{d}z \, e^{(s-1)L(t)} G_{\rho}(0, z) \\ &= \frac{i\Gamma(1-s)}{4 \sin(8\pi s)} \int_{\gamma(t)} \mathrm{d}z \, e^{sL(t)}\\ &= \frac{i\Gamma(1-s)}{4 \sin(8\pi s)} \int_{\gamma(t)} \mathrm{d}z \, [e^{sL(t)} - 1] \\ &\to \frac{i}{32\pi} \int_{\gamma(t)} \mathrm{d}z \, L(t) \\ &= \frac{i}{32\pi} \left[ \int_{-1}^{1} [\log(1-z^2) - i\pi] \, \mathrm{d}z - \int_{-1}^{1} [\log(1-z^2) + i\pi] \, \mathrm{d}z \right] \\ &= \frac{1}{8}. \tag{2.5}\label{e:2.5} \end{align*} $$ Similarly, the remaining terms for $\eqref{e:2.4}$ can be recast as $$ \begin{align*} &- \frac{\Gamma(1-s)}{4\pi \sin(8\pi s)} \int_{\gamma(t)} \mathrm{d}z \, e^{(s-1)L(t)} [G_{\rho}(s, z) - G_{\rho}(1, z) ] \\ &\to - \frac{1}{32\pi^2} \int_{\gamma(t)} \mathrm{d}z \, e^{-L(t)} \frac{\partial G_{\rho}}{\partial s} (s, z) \\ &= \frac{1}{32\pi^2} \int_{\gamma(t)} \mathrm{d}z \, \frac{1}{1 - z^2} \int_{C(\rho)} \mathrm{d}y \, \frac{\sqrt{1 + 2zy + y^2}}{y^3} \\ &\hspace{5em} \times \left[ 8\log y - \log(1-y^2) - \frac{5}{2}\log(1+2zy+y^2) \right] \tag{2.6}\label{e:2.6} \end{align*} $$ We can first remove the integral with respect to $z$ by invoking the residue theorem. Indeed,

Observation. For any function $f(z)$ analytic near $[-1, 1]$, $$\begin{align*} \int_{\gamma(t)} \mathrm{d}z \, \frac{f(z)}{1 - z^2} = - i\pi [f(1) + f(-1)]. \end{align*}$$

From this, we get $$ \begin{align*} \eqref{e:2.6} &= -\frac{i}{32\pi} \int_{C(\rho)} \mathrm{d}y \, \frac{1 + y}{y^3} \left( 8\log y - \log(1-y^2) - 5\log(1+y) \right) \\ &\qquad -\frac{i}{32\pi} \int_{C(\rho)} \mathrm{d}y \, \frac{1-y}{y^3} \left[ 8\log y - \log(1-y^2) - 5\log(1-y) \right] \\ &\to -\frac{i}{32\pi} (-22i\pi) = -\frac{11}{16} \qquad \text{as} \quad \rho \to 1^-. \end{align*} $$ Consequently, it follows that $$ B_1(0) = \eqref{e:2.5} + \eqref{e:2.6} = -\frac{9}{16}. $$ Hence, $\eqref{e:2.1}$ shows that the contribution from $\eqref{e:1b}$ is given by $$ \begin{align*} \lim_{s\to 0} \int_{0}^{\infty} \mathrm{d}y \, \frac{\eqref{e:1b}}{y^s(y+1)^s(y+2)^s} &= \frac{2^{-2}\Gamma(-\frac{1}{2})}{\Gamma(\frac{1}{2})} B_1(0) = \frac{9}{32}. \tag{2.7}\label{e:2.7} \end{align*} $$

3. Contribution from $\eqref{e:1a}$

Now we tackle $\eqref{e:1a}$: $$ \begin{align*} &\int_{0}^{\infty} \mathrm{d}y \, \frac{\eqref{e:1a}}{y^s(y+1)^s(y+2)^s} \\ &= - \frac{\sin(\pi s)}{\sin(2\pi s)} \int_{0}^{1} \frac{\mathrm{d}x}{x^s (1-x)^s} \int_{0}^{\infty} \frac{\mathrm{d}y}{y^s(y+1)^s(y+2)^s (y + 2x)^s (y+x)^s (y+1+x)^s} \end{align*} $$ As claimed at the beginning, this is much easier. Indeed, noting that the leftmost zero of the denominator is $y = -2$, for $\Re(s) > \frac{1}{6}$ we have $$ \int_{\infty}^{(-2)^+} \frac{\mathrm{d}y}{(-y)^s(-y-1)^s(-y-2)^s (-y- 2x)^s (-y-x)^s (-y-1-x)^s} = 0, $$ where the integral is taken along the Hankel contour begins at $\infty$, circles the point $-2$ once in the positive direction, and returns to $\infty$. Simplifying, $$ \begin{align*} & \int_{0}^{\infty} \frac{\mathrm{d}y}{y^s(y+1)^s(y+2)^s (y + 2x)^s (y+x)^s (y+1+x)^s} \\ &= -\frac{1}{\sin(6\pi s)} \int_{0}^{2} \frac{\sin(\#[\{1, 2, x, x+1, 2x\} \cap [y, \infty)]\pi s) \, \mathrm{d}y}{y^s |(1-y) (2-y) (x-y) (x+1-y) (2x-y)|^s}, \end{align*} $$ where $\#A$ denotes the number of elements of the set $A$. Note that this function is meromorphic for $\Re(s) < 1$. Plugging this back, it follows that $\eqref{e:1a}$ defines a meromorphic function for $\Re(s) < 1$. Now letting $s \to 0$, $$ \begin{align*} &\lim_{s\to 0}\int_{0}^{\infty} \mathrm{d}y \, \frac{\eqref{e:1a}}{y^s(y+1)^s(y+2)^s} \\ &= \frac{1}{12} \int_{0}^{1} \mathrm{d}x \int_{0}^{2} \mathrm{d}y \, \#[\{1, 2, x, x+1, 2x\} \cap [y, \infty)] \\ &= \frac{1}{12} \int_{0}^{1} \mathrm{d}x \int_{0}^{2} \mathrm{d}y \, \bigl( \mathbf{1}[1 > y] + \mathbf{1}[2 > y] + \mathbf{1}[x > y] + \mathbf{1}[x+1 > y] + \mathbf{1}[2x > y] \bigr) \\ &= \frac{1}{12} \int_{0}^{1} \mathrm{d}x \, (1 + 2 + x + (x+1) + 2x) \\ &= \frac{1}{2}. \tag{3.1} \label{e:3.1} \end{align*} $$

4. Conclusion

In summary, the function $$ F(s) = \int_{0}^{\infty}\int_{0}^{\infty} \frac{\mathrm{d}x\mathrm{d}y}{x^s (x+1)^s y^s (y+1)^s (y+2)^s (2x + y + 2)^s (x+y+1)^s (x+y+2)^s} $$ extends to a meromorphic function at least for $\Re(s) < \frac{1}{3}$ and satisfies $$ F(0) = \eqref{e:2.7} + \eqref{e:3.1} = \frac{25}{32}. $$

Answer 2

I think your claim is true: the analytic continuation of $F$ is defined at $\mathbb Z^- := \{0, -1, -2, \dots\}$ and takes on rational values there. In this answer I will show how to construct this analytic continuation, and I will give a (partially rigorous) method to compute $F(0) = 25/32$. I don't think there is a simple formula for $F(-n)$, but there is a somewhat practical algorithm to compute it, which could be implemented in Sage or Mathematica.

NB: I edited my answer, improving the exposition and the proof method, fixing computational errors, and sketching a partial proof of the claim. There are a couple unproven assertions in bold, and more importantly, there is a heuristic I do not know how to prove or make rigorous. Any suggestions are welcome.

Let $f(x,y,s)$ be your integrand. Define:

\begin{align} R_{00} &:= [0,1]\times [0,1] \\ R_{10} &:= [1,+\infty)\times [0,1] \\ R_{01} &:= [0,1]\times [1,+\infty) \\ R_{11} &:= [1,+\infty)\times [1,+\infty) \\ R_{11}^+ &:= R_{11} \cap \{x\ge y\} \\ R_{11}^- &:= R_{11} \cap \{x\le y\} \\ \end{align} Then, define $$I_{00}(s) := \iint_{R_{00}} f(x,y,s)\,dx\,dy$$ and so on.

Now, clearly $F(s) = I_{00}(s) + I_{10}(s) + I_{01}(s) + I_{11}^+(s) + I_{11}^-(s)$. For each summand we will:

• analytically continue it;
• prove that the continuation is rational at each non-positive integer;
• compute its value at $s = 0$.

---

The integral $I_{00}(s)$ is convergent for $\Re(s) < 1$, and the integrand is a polynomial for $s\in \mathbb Z^-$, so the integral is rational there. Moreover, $I_{00}(0) = \boxed{1}$.

---

The integral $I_{10}(s)$ is convergent for $1/5 < \Re(s) < 1$, and will serve as the introduction to the main technique.

Notice that in $R_{10}$, we can expand $f$ for $x\to +\infty$: \begin{align} f(x,y,s) &= (2x^5y(y+1)(y+2))^{-s}\left[\left(1 + \frac{1}{x}\right)^{-s} \left(1 + \frac{y+1}{x}\right)^{-s} \left(1 + \frac{y+2}{2x}\right)^{-s} \left(1 + \frac{y+2}{x}\right)^{-s}\right] \\ &= (2x^5y(y+1)(y+2))^{-s} \left[ \sum_{k = 0}^n a_k(y, s)x^{-k}+ O(x^{-n-1})\right], \end{align} where, crucially, the $a_k(y, s)$ are polynomials in $y, s$ with rational coefficients. Moreover, if $s \in \mathbb Z^-$, the series is a product of four finite binomial expansions, so $a_k(y,s) = 0$ for $s \in \{0, -1, \dots, -\lfloor(k-1)/4\rfloor\}$.

Let us now integrate the truncated series by first assuming $1/5 < \Re(s) < 1$:

\begin{align} &\iint_{R_{10}} (2x^5y(y+1)(y+2))^{-s} \left[ \sum_{k = 0}^n a_k(y, s)x^{-k}\right] \, dx\, dy \\ &= \sum_{k = 0}^n\int_{1}^{+\infty} x^{-5s-k}\, dx \int_0^1 (2y(y+1)(y+2))^{-s}a_k(y, s) \, dy \\ &= \sum_{k = 0}^n \frac{1}{5s+k-1} \int_0^1 (2y(y+1)(y+2))^{-s}a_k(y, s) \, dy =: A_n(s). \end{align} The expression for $A_n(s)$ defines a meromorphic function on the set $\{\Re(s)<1\}$, with the property that $A_n(s) \in \mathbb Q$ for all $s \in \mathbb Z^-$. Indeed, any simple pole at $s = (1-k)/5 \in \mathbb Z^-$ is cancelled out by the factor $a_k(y,s)$. Hence, we can write $$ I_{10}(s) = A_n(s) + \iint_{R_{10}}\left[f(x,y,s) - (2x^5y(y+1)(y+2))^{-s} \left[ \sum_{k = 0}^n a_k(y, s)x^{-k}\right]\right] \, dx\, dy. $$ Note that now the integrand is $O(x^{-n-5s-1})$, so $I_{10}(s)$ is now defined and meromorphic on $\{-n/5 < \Re(s) < 1\}$. Moreover, at $s = 0, -1, \dots, -\lfloor n/4 \rfloor$, the integrand vanishes and $A_n(s)$ is rational. By choosing $n$ large enough, we see that $I_{10}(s)$ is rational for all $s \in \mathbb Z^-$.

To compute $I_{10}(0)$, it suffices to take $n=1$. We have \begin{align} a_0(y,s) &= 1, \\ a_1(y,s) &= -\frac{5}{2}s(y+2). \end{align} Hence, \begin{align} I_{10}(0) = A_1(0) &= -1\int_{0}^{1}a_0(y,s)\,dy + \frac{1}{5}\int_{0}^{1}[a_1(y,s)/s]\big|_{s=0}\,dy \\ &= -1 - \frac{1}{5}\int_{0}^{1}\frac{5}{2}(y+2)\,dy \\ &= -1-5/4 = \boxed{-9/4}. \end{align}

---

The integral $I_{01}(s)$ can be analytically continued to a meromorphic function on $\{\Re(s) < 1\}$ in a very similar way. First expand $$ f(x,y,s) = (y^6x(x+1))^{-s}\left[ \sum_{k = 0}^n b_k(x, s)y^{-k}+ O(y^{-n-1})\right]. $$ The truncated series integrates to $$ B_n(s) := \sum_{k=0}^n \frac{1}{6s+k-1} \int_0^1 (x(x+1))^{-s}b_k(x, s)\,dx, $$ which is meromorphic for $\Re(s) < 1$ and rational for $s \in \mathbb Z^-$. This implies that $I_{01}(s)$ is rational for $s \in \mathbb Z^-$.

As before, we take $n=1$ and compute \begin{align} b_0(x,s) &= 1, \\ b_1(x,s) &= -4s(x+2), \end{align} \begin{align} I_{01}(0) = B_1(0) &= -1\int_{0}^{1}b_0(x,s)\,dx + \frac{1}{6}\int_{0}^{1}[b_1(x,s)/s]\big|_{s=0}\,dx \\ &= -1 - \frac{1}{6}\int_{0}^{1}4(x+2)\,dx \\ &= -1-5/3 = \boxed{-8/3}. \end{align}

---

The integral $I_{11}^+(s)$ is similar, but first we have to get rid of the annoying factors $(x+y+2)$ etc. by substituting $x=yt$:

\begin{align} I_{11}^+(s) &= \int_1^{+\infty}dy \int_y^{+\infty}dx\,(xy(x+1)(y+1)(y+2)(x+y+1)(2x+y+2)(x+y+2))^{-s} \\ &= \int_1^{+\infty}dy \int_1^{+\infty}y\,dt\,(y^2t(yt+1)(y+1)(y+2)(yt+y+1)(2yt+y+2)(yt+y+2))^{-s}. \end{align} Let \begin{align} R(t,y)&:= \left(1+\frac{1}{yt}\right) \left(1+\frac{y+1}{yt}\right) \left(1+\frac{y+2}{2yt}\right) \left(1+\frac{y+2}{yt}\right), \\ S(y)&:=\left(1+\frac{1}{y}\right)\left(1+\frac{2}{y}\right). \end{align} Expand \begin{align} R(t,y)^{-s} &:= \sum_{k=0}^n c_k^+(y,s)t^{-k} + O(t^{-n-1}), \\ c_k^+(y,s) S(y)^{-s} &:= \sum_{j=0}^m c_{k,j}^+(s)y^{-j} + O(y^{-m-1}), \end{align} and notice that if $s \in \mathbb Z^-$, then $R(t,y)$ is a finite sum, so $c_k^+(y,s) = 0$ for $s = 0, -1, \dots, -\lfloor (k-1)/4\rfloor$. Moreover, $c_k^+(y,s) \in y^{-k}\mathbb Q[y,s]$, with $\deg_y (y^kc_k^+(y,s)) = k$, since the coefficient of $(-s)^k$ in $y^kc_k^+(y,s)$ is a polynomial in $y$ with positive coefficients. Hence, if $s \in \mathbb Z^-$, then $c_k^+(y,s) S(y)^{-s}$ involves only the terms $1, y^{-1}, \dots, y^{-k+2s}$, and $c^+_{k,j}(s) = 0$ for $s=0, -1, \dots, -\max(\lfloor (k-1)/4 \rfloor,\lfloor (j-1)/6 \rfloor)$.

We have \begin{align} I_{11}^+(s) &= \iint_{R_{11}} y(2y^8t^5)^{-s}S(y)^{-s}R(t,y)^{-s}\,dy\,dt \\ &= \iint_{R_{11}} y(2y^8t^5)^{-s}S(y)^{-s}\left[R(t,y)^{-s}-\sum_{k=0}^n c_k^+(y,s)t^{-k}\right]\,dy\,dt \\ &+ \sum_{k=0}^n\iint_{R_{11}} y(2y^8t^5)^{-s}t^{-k}\left[S(y)^{-s} c_k^+(y,s)- \sum_{j=0}^m c_{k,j}^+(s)y^{-j} \right] \,dy\,dt \\ &+ \sum_{k=0}^n\sum_{j=0}^m \iint_{R_{11}} y(2y^8t^5)^{-s}t^{-k}y^{-j}c_{k,j}^+(s) \,dy\,dt \\ &= \iint_{R_{11}} y(2y^8t^5)^{-s}S(y)^{-s}\left[R(t,y)^{-s}-\sum_{k=0}^n c_k^+(y,s)t^{-k}\right]\,dy\,dt \tag{$1$} \\ &+ \sum_{k=0}^n\frac{2^{-s}}{5s+k-1}\int_{1}^{+\infty} y^{1-8s}\left[S(y)^{-s} c_k^+(y,s)- \sum_{j=0}^m c_{k,j}^+(s)y^{-j} \right] \,dy \tag{$2$} \\ &+ \sum_{k=0}^n\sum_{j=0}^m \frac{2^{-s}c_{k,j}^+(s)}{(5s+k-1)(8s+j-2)} \tag{$3$} \end{align}

The integral $(3)$ is a meromorphic function, which takes on finite rational values on $\mathbb Z^-$, because any poles are canceled out by $c_{k,j}^+(s)$ [I have not checked this, but it should be true. It involves looking at the binomial coefficients in the series expansion.] Moreover, for any given $s \in \mathbb Z^-$, the integral $(2)$ vanishes for $m$ large enough.

The following is not correct:

Hence, in the case $s = 0$, we have $$ I_{11}^+(0) := C_{1,2}^+(0) = \left[\sum_{k=0}^1 \sum_{j=0}^2 \frac{c_{k,j}^+(s)}{(5s+k-1)(8s+j-2)}\right]_{s = 0}. $$ Using Mathematica, we find \begin{align} c_{0,0}^+(s) &= 1 & c_{1,0}^+(s) &= -\frac{5s}{2} \\ c_{0,1}^+(s) &= -3s & c_{1,1}^+(s) &= -\frac{5s}{2}(2-3s) \\ c_{0,2}^+(s) &= \frac{s}{2}(5+9s) & c_{1,2}^+(s) &= -\frac{5s^2}{4}(9s-7) \\ \end{align} whence $I_{11}^+(0) = \boxed{53/32}$.

What actually happens (as computed by Mathematica) is that $$ C_{n,2}^+(0) = \left(\frac{3}{16}, \frac{53}{32}, \frac{67}{32}, \frac{7}{4}, \frac{133}{64}, \frac{225}{128}, \frac{531}{256}, \frac{901}{512}, \frac{2123}{1024}, \frac{3605}{2048}, \dots \right), $$ $$ C_{n,2}^+(0) = \begin{cases} 3/16 & n = 0 \\ 53/32 & n = 1 \\ \frac{1}{12}\left(23 + \frac{15}{8} (-1)^n + \left(-\frac{1}{2}\right)^n\right) & n \ge 2 \end{cases} $$ and hence $\frac{1}{6}(C_{n,2}^+(0) + 3C_{n+1,2}^+(0) + 2C_{n+1,2}^+(0)) \equiv \boxed{23/12}$ for $n\ge 2$. I do not know how to prove this or even how to make this rigorous, but I suspect that the coefficients $(1,3,2)$ have to do with those of $S(y)$ or $R(t,y)\big|_{t\to +\infty}$.

---

A similar procedure "computes" $I_{11}^-(0)$, starting with a substitution $y = xt$, and ending with a "limit regularization". The inferred value is $I_{11}^-(0) = \boxed{89/32}$.

---

To summarize, if we can fix the "limit regularization" steps, we have that $F(s)$ is well defined, analytic and rational at every $s \in \mathbb Z^-$, and moreover $$ F(0) = 1 - 9/4 - 8/3 + 23/12 + 89/32 = \boxed{25/32}. $$ This is in agreement with your guess, which indicates that the general method can probably be made rigorous.