Is this Euler-Mascheroni constant calculation from double integrals a true identity?

Question

A prime number is a number that is only divisible by itself and one, that is the number of divisors of a prime number is equal to $2$. One way to illustrate this is to plot a matrix such that if the column index (1,2,3,...) divides the row index (1,2,3,...) then a black square is drawn at that column and row index intersection, like this:

As a matrix this has the definition in Mathematica:

MatrixForm[Table[Table[If[Mod[n, k] == 0, 1, 0 ], {k, 1, 12}], {n, 1, 12}]]

which gives the table below:

$$\small \left( \begin{array}{cccccccccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 1 & 1 & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \end{array} \right)$$

The number of divisors of n, also called tau, is the row sums of this matrix and starts: 1, 2, 2, 3, 2, 4, 2, 4, 3, 4, 2, 6 ...

So at a row equal to a prime number there are only two black dots and the row sums in the above matrix are equal to 2.

To find the ordinary generating function of the number of divisors of n we consider table $a$ times $b$, or row index times column index:

MatrixForm[Table[Table[a*b, {b, 1, 12}], {a, 1, 12}]]

$$\small \left( \begin{array}{cccccccccccc} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 \\ 2 & 4 & 6 & 8 & 10 & 12 & 14 & 16 & 18 & 20 & 22 & 24 \\ 3 & 6 & 9 & 12 & 15 & 18 & 21 & 24 & 27 & 30 & 33 & 36 \\ 4 & 8 & 12 & 16 & 20 & 24 & 28 & 32 & 36 & 40 & 44 & 48 \\ 5 & 10 & 15 & 20 & 25 & 30 & 35 & 40 & 45 & 50 & 55 & 60 \\ 6 & 12 & 18 & 24 & 30 & 36 & 42 & 48 & 54 & 60 & 66 & 72 \\ 7 & 14 & 21 & 28 & 35 & 42 & 49 & 56 & 63 & 70 & 77 & 84 \\ 8 & 16 & 24 & 32 & 40 & 48 & 56 & 64 & 72 & 80 & 88 & 96 \\ 9 & 18 & 27 & 36 & 45 & 54 & 63 & 72 & 81 & 90 & 99 & 108 \\ 10 & 20 & 30 & 40 & 50 & 60 & 70 & 80 & 90 & 100 & 110 & 120 \\ 11 & 22 & 33 & 44 & 55 & 66 & 77 & 88 & 99 & 110 & 121 & 132 \\ 12 & 24 & 36 & 48 & 60 & 72 & 84 & 96 & 108 & 120 & 132 & 144 \end{array} \right)$$

By stretching out this table in the downward direction we get this table:

MatrixForm[Table[Table[If[Mod[a, b] == 0, a, 0], {b, 1, 12}], {a, 1, 12}]]

$$\small \left( \begin{array}{cccccccccccc} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 2 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 3 & 0 & 3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 4 & 4 & 0 & 4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 5 & 0 & 0 & 0 & 5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 6 & 6 & 6 & 0 & 0 & 6 & 0 & 0 & 0 & 0 & 0 & 0 \\ 7 & 0 & 0 & 0 & 0 & 0 & 7 & 0 & 0 & 0 & 0 & 0 \\ 8 & 8 & 0 & 8 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 \\ 9 & 0 & 9 & 0 & 0 & 0 & 0 & 0 & 9 & 0 & 0 & 0 \\ 10 & 10 & 0 & 0 & 10 & 0 & 0 & 0 & 0 & 10 & 0 & 0 \\ 11 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 11 & 0 \\ 12 & 12 & 12 & 12 & 0 & 12 & 0 & 0 & 0 & 0 & 0 & 12 \end{array} \right)$$

This looks very much like the divisorplot at the beginning. To find the generating function for the number of divisors of $n$ we then introduce the variable $x$ like this:

Clear[x]
MatrixForm[
 Table[Table[If[Mod[a, b] == 0, x^a, 0], {b, 1, 12}], {a, 1, 12}]]

$$\small \left( \begin{array}{cccccccccccc} x & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^2 & x^2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^3 & 0 & x^3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^4 & x^4 & 0 & x^4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^5 & 0 & 0 & 0 & x^5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^6 & x^6 & x^6 & 0 & 0 & x^6 & 0 & 0 & 0 & 0 & 0 & 0 \\ x^7 & 0 & 0 & 0 & 0 & 0 & x^7 & 0 & 0 & 0 & 0 & 0 \\ x^8 & x^8 & 0 & x^8 & 0 & 0 & 0 & x^8 & 0 & 0 & 0 & 0 \\ x^9 & 0 & x^9 & 0 & 0 & 0 & 0 & 0 & x^9 & 0 & 0 & 0 \\ x^{10} & x^{10} & 0 & 0 & x^{10} & 0 & 0 & 0 & 0 & x^{10} & 0 & 0 \\ x^{11} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & x^{11} & 0 \\ x^{12} & x^{12} & x^{12} & x^{12} & 0 & x^{12} & 0 & 0 & 0 & 0 & 0 & x^{12} \end{array} \right)$$

The generating function for the number of divisors of $n$ is then the infinite double sum which is the sum of the infinite matrix above:

$$\sum\limits_{n=1}^{\infty} d(n)x^{n} = \sum\limits_{a=1}^{\infty} \sum\limits_{b=1}^{\infty} x^{a b}$$ $$= 1\, x^{1} + 2\, x^{2} + 2\, x^{3} + 3\, x^{4} + 2\, x^{5} + 4\, x^{6} + 2\, x^{7} + 4\, x^{8} + 3\, x^{9} + 4\, x^{10} + 2\, x^{11} + 6\, x^{12} + ...$$

where $d(n)$ = the number of divisors of $n$

A sum $\sum$ is a sum of distances - usually under a function, while an integral $\int$ is a sum of areas under a function, and a double integral $\int \int$ is a sum of volumes and so on... Often in mathematics one wants to compare a sum to an integral, or at least so I have heard. In this case we have a double sum:

$$\sum\limits_{a=1}^{\infty} \sum\limits_{b=1}^{\infty} x^{a b}$$

Considering the double integral instead we have:

$$\int\limits_{a=1}^{\infty} \int\limits_{b=1}^{\infty} x^{a b}$$

As suggested by Chris's wise sister in the chat room about two weeks ago, the similar triple integral evaluates faster at $x=\frac{1}{e}$:

According to Mathematica:

Integrate[Integrate[1/Exp[1]^(a*b), {a, 1, Infinity}], {b, 1, Infinity}]

the double integral then evaluates to:

$$\int\limits_{a=1}^{\infty} \int\limits_{b=1}^{\infty} \frac{1}{e^{a b}}=-\text{Ei}(-1)$$

Or as Mathematica calls it, the exponential integral. Reading about the exponential integral in Wikipedia that passing a logarithm as argument to it one gets the logarithmic integral and remembering that the logarithmic integral has to do with the prime number theorem, I did some experimenting in Mathematica with the program:

Monitor[MatrixForm[
  Table[Table[
    Integrate[
      Integrate[Exp[a*b], {a, k, Infinity}, 
       GenerateConditions -> False], {b, 0, n}, 
      GenerateConditions -> False] - 
     Integrate[
      Integrate[Exp[a*b], {a, 0, n}, GenerateConditions -> False], {b,
        k, Infinity}, GenerateConditions -> False], {n, 1, 4}], {k, 0,
     4}]], k]

and came to the conclusion that for $n \gt 0$ and $k \geq 0$, this appears to be an identity:

$$\gamma = \int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db$$

where $\gamma$ = Euler-Mascheroni constant, approximately equal to 0.57721566490153286060651209...

My question is: Why is this identity true?

---

I also noted that:

For $k \lt 0$ and $n \lt 0$, it appears to be true that:

$$\gamma = \Re(\int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db)$$

$$-\pi = \Im(\int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db)$$

But this is probably less interesting.

---

As a starting point I believe it has something to do with the fact:

-Integrate[Integrate[Exp[a*b], {a, 0, 1}, GenerateConditions -> False], {b, 0, 
   Infinity}, GenerateConditions -> False]

$$\gamma = -\int_0^{\infty } \left(\int_0^1 e^{a b} \, da\right) \, db$$

and that it somehow can be separated out from both of the double integrals: $$\int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db$$ $$\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db$$

depending on the values of $n$ and $k$ since:

Monitor[MatrixForm[
  Table[Table[
    Integrate[
      Integrate[Exp[a*b], {a, k, Infinity}, 
       GenerateConditions -> False], {b, 0, n}, 
      GenerateConditions -> False] - 
     0*Integrate[
       Integrate[Exp[a*b], {a, 0, n}, 
        GenerateConditions -> False], {b, k, Infinity}, 
       GenerateConditions -> False], {n, 1, 4}], {k, 0, 4}]], k]
Monitor[MatrixForm[
  Table[Table[
    0*Integrate[
       Integrate[Exp[a*b], {a, k, Infinity}, 
        GenerateConditions -> False], {b, 0, n}, 
       GenerateConditions -> False] - 
     Integrate[
      Integrate[Exp[a*b], {a, 0, n}, GenerateConditions -> False], {b,
        k, Infinity}, GenerateConditions -> False], {n, 1, 4}], {k, 0,
     4}]], k] 

for $n \geq 1$ and $k \geq 0$ we have for the first double integral:

$$\int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db$$

(for $n=1,2,3,4...$ and $k=0,1,2,3,4...$):

$$\small \left( \begin{array}{cccc} 0 & -\log (2) & -\log (3) & -\log (4) \\ \gamma -\text{Ei}(1) & \gamma -\text{Ei}(2) & \gamma -\text{Ei}(3) & \gamma -\text{Ei}(4) \\ -\text{Ei}(2)+\gamma +\log (2) & -\text{Ei}(4)+\gamma +\log (2) & -\text{Ei}(6)+\gamma +\log (2) & -\text{Ei}(8)+\gamma +\log (2) \\ -\text{Ei}(3)+\gamma +\log (3) & -\text{Ei}(6)+\gamma +\log (3) & -\text{Ei}(9)+\gamma +\log (3) & -\text{Ei}(12)+\gamma +\log (3) \\ -\text{Ei}(4)+\gamma +\log (4) & -\text{Ei}(8)+\gamma +\log (4) & -\text{Ei}(12)+\gamma +\log (4) & -\text{Ei}(16)+\gamma +\log (4) \end{array} \right)$$

and for the second double integral:

$$-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db)$$

(for $n=1,2,3,4...$ and $k=0,1,2,3,4...$):

$$\small \left( \begin{array}{cccc} \gamma & \gamma +\log (2) & \gamma +\log (3) & \gamma +\log (4) \\ \text{Ei}(1) & \text{Ei}(2) & \text{Ei}(3) & \text{Ei}(4) \\ \text{Ei}(2)-\log (2) & \text{Ei}(4)-\log (2) & \text{Ei}(6)-\log (2) & \text{Ei}(8)-\log (2) \\ \text{Ei}(3)-\log (3) & \text{Ei}(6)-\log (3) & \text{Ei}(9)-\log (3) & \text{Ei}(12)-\log (3) \\ \text{Ei}(4)-\log (4) & \text{Ei}(8)-\log (4) & \text{Ei}(12)-\log (4) & \text{Ei}(16)-\log (4) \end{array} \right)$$

I was not entirely clear in the use of integration indices $n$ and $k$. Usually $n$ is the row index and $k$ is the column index, while in the matrices above they have been permuted.

---

Edit 10.6.2013: I should have used the 'PrincipalValue' command instead of the 'GenerateConditions' command, like this:

t = 1;
Table[Integrate[
  Integrate[Exp[a*b] + 1/b/t, {a, 0, t}], {b, -Infinity, n}, 
  PrincipalValue -> True], {n, 1, 4}]
Table[Integrate[
  Integrate[Exp[a*b] + 1/b/t, {a, 0, t}], {b, -Infinity, Log[n]}, 
  PrincipalValue -> True], {n, 1, 4}]

$$\text{Ei}(n \cdot t) = \int_{-\infty }^n \left(\int_0^t \left(e^{a b}+\frac{1}{b \cdot t}\right) \, da\right) \, db$$

And the relation from Wikipedia I wanted explore is: $$\text{Ei}(\log (n)) = \text{li}(n)$$

Edit 11.6.2013:

This allows us to write:

X=5;
Monitor[Table[
  Integrate[
   Integrate[Exp[a*b] + n/b, {a, 0, 1/n}], {b, -Infinity, Log[X]}, 
   PrincipalValue -> True], {n, 1, 6}], n]

$$\text{li}(X^{\frac{1}{n}}) = \int_{-\infty }^{\text{Log}[X]} \left(\int_0^{\frac{1}{n}} \left(e^{a b}+\frac{n}{b}\right) \, da\right) \, db$$

Which simplifies to:

X = 5;
Monitor[Table[
  Integrate[E^(b/n)/b, {b, -Infinity, Log[X]}, 
   PrincipalValue -> True], {n, 1, 6}], n]

$$\text{li}(X^{\frac{1}{n}}) = \int_{-\infty }^{\text{Log}[X]} \frac{e^{b/n}}{b} \, db$$

Answer 1

What justifies diligently the step $\sum\limits_{a=1}^{\infty} \sum\limits_{b=1}^{\infty} x^{a b} \rightarrow \int\limits_{a=1}^{\infty} \int\limits_{b=1}^{\infty} x^{a b}$; this seems likely very vague "what is valid for distances ... valid for areas ... valid for volumes..." so might one argument and extend to higher dimmensions $\int\int\int\int...$ So far cannot follow to justify identity.

Continue that there is no reason to me to believe that this identity works, because if then according to Hadjicostas/Sondow you would have:

$$ \int_0^n \left(\int_k^{\infty } e^{a b} \, da\right) \, db-\int_k^{\infty } \left(\int_0^n e^{a b} \, da\right) \, db=-\int^1_0\int^1_0\frac{1-a}{(1-ab)\log (ab)} da\; db = \gamma$$

isn't it?