I am seeking a closed-form solution for this double sum:
\begin{eqnarray*} \sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(\color{blue}{3}n+m)}= ?. \end{eqnarray*}
I will turn it into $3$ tough integrals in a moment. But first I will state some similar results:
\begin{eqnarray*} \sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(n+m)} &=& 2 \zeta(3) \\ \sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(\color{blue}{2}n+m)} &=& \frac{11}{8} \zeta(3) \\ \sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \frac{1}{nm(\color{blue}{4}n+m)} &=& \frac{67}{32} \zeta(3) -\frac{G \pi}{2}. \\ \end{eqnarray*}
where $G$ is the Catalan constant. The last result took some effort ...
Now I know most of you folks prefer integrals to sums, so lets turn this into an integral. Using
\begin{eqnarray*} \frac{1}{n} &=& \int_0^1 x^{n-1} dx\\ \frac{1}{m} &=& \int_0^1 y^{m-1} dy\\ \frac{1}{3n+m} &=& \int_0^1 z^{3n+m-1} dz \\ \end{eqnarray*} and summing the geometric series, we have the following triple integral \begin{eqnarray*} \int_0^1 \int_0^1 \int_0^1 \frac{z^3 dx dy dz}{(1-xz^3)(1-yz)}. \end{eqnarray*}
Now doing the $x$ and $y$ integrations we have \begin{eqnarray*} I=\int_0^1 \frac{\ln(1-z) \ln(1-z^3)}{z} dz. \end{eqnarray*}
Factorize the argument of the second logarithm ...
\begin{eqnarray*} I= \underbrace{\int_0^1 \frac{\ln(1-z) \ln(1-z)}{z} dz}_{=2\zeta(3)} + \int_0^1 \frac{\ln(1-z) \ln(1+z+z^2)}{z} dz. \end{eqnarray*}
So if you prefer my question is ... find a closed form for:
\begin{eqnarray*} I_1 = - \int_0^1 \frac{\ln(1-z) \ln(1+z+z^2)}{z} dz. \end{eqnarray*}
Integrating by parts gives:
\begin{eqnarray*} I_1 = - \int_0^1 \frac{\ln(z) \ln(1+z+z^2)}{1-z} dz + \int_0^1 \frac{(1+2z)\ln(z) \ln(1-z)}{1+z+z^2} dz. \end{eqnarray*}
and let us call these integrals $I_2$ and $I_3$ respectively.
All $3$ of these integrals are not easy for me to evaluate and any help with their resolution will be gratefully received.
Let's introduce a function:
$$S(x,y)=\sum_{n=1}^\infty \sum_{m=1}^\infty \frac{x^n y^m}{ n m (3n +m)}$$
Assuming $|x|<1$ and $|y|<1$ we avoid any convergence problems, and can use partial fractions:
$$S(x,y)=\sum_{n=1}^\infty \sum_{m=1}^\infty \frac{x^n y^m}{ n m^2}-\sum_{n=1}^\infty \sum_{m=1}^\infty \frac{x^n y^m}{ m^2(n+ \frac{1}{3} m)}$$
$$S(x,y)=-\log(1-x) \text{Li}_2(y) -x \sum_{m=1}^\infty \frac{y^m}{ m^2} \Phi \left(x,1,\frac13 m+1 \right)$$
Let's use the integral representation of the Lerch transcendent:
$$\Phi \left(x,1,\frac13 m+1 \right)= \int_0^\infty \frac{e^{-(1+\frac13 m)t} ~dt}{1-x e^{-t}}$$
Summation under the integral gives us:
$$S(x,y)=-\log(1-x) \text{Li}_2(y) -x \int_0^\infty \text{Li}_2 \left(y e^{-t/3} \right) \frac{dt}{e^t-x}$$
So we can assume:
Which does seem to work numerically, though of course is quite hard to evaluate symbolically.
Numerical check:
Compare to exact expression:
It may be that $x=y$ is not the best choice for the limit. For example, we can assume $x=y^a$ where $a$ is some real number. A good choice may lead to a better numerical convergence or even a closed form.
Using the dilogarithm properties, we have:
$$F(x,y)=-x \int_0^\infty \text{Li}_2 \left(y e^{-t/3} \right) \frac{dt}{e^t-x}=x \int_0^\infty \int_0^1 \frac{\log(1-e^{-t/3} y u) du dt}{u (e^t-x)}$$
Let's change the variable:
$$e^{-t}=v \\ t=- \log v$$
$$F(x,y)=x \int_0^1 \int_0^1 \frac{\log(1- y u v^{1/3}) du dv}{u (1-x v)}$$
Let's take:
$$y=x^{1/3}$$
We have:
$$F(x,y)=x \int_0^1 \int_0^1 \frac{\log(1- u (xv)^{1/3}) du dv}{u (1-x v)}$$
$$v=w/x$$
$$F(x,y)=\int_0^x \int_0^1 \frac{\log(1- u w^{1/3}) du dw}{u (1-w)}$$
$$F(x,y)=-\int_0^x \frac{\text{Li}_2 (w^{1/3}) dw}{1-w}$$
So, there's a neater expression for the limit:
This allows a simple generalization:
$$\sum_{n=1}^\infty \sum_{m=1}^\infty \frac{x^n y^m}{ n m (an +m)}=\lim_{x \to 1} \left[-\log(1-x) \text{Li}_2(x^{1/a}) -\int_0^x \frac{\text{Li}_2 (w^{1/a}) dw}{1-w} \right]$$
Which checks out numerically with the examples from the OP.
I wonder if we can somehow use L'Hospital here to deal with the integral and get a closed form for the limit.
Integration by parts could also work.