An Euler-Mascheroni double sum An interesting  representation of the Euler-Mascheroni constant
$$ \gamma~=~ \lim \limits_{n\to \infty} \sum \limits_{k,s=1}^n \frac{s-k}{k\left( s\,n +k\right)},\label{1}\tag{$*$}$$
can be proved using the relatively recent identity of Macys (Mat. Zametki, 94/5, 695–701, 2013), but it also can be reformulated  to
almost mimic a Riemann sum over the unit square $[0,1]^2$
$$ \gamma~=~ \lim \limits_{n\to \infty} \sum \limits_{x,y=1/n\,(step: 1/n)}^1 \frac{x-y}{y\,\left( x +y/n\right)} \,\left(\frac1n\right)^2,$$
and there are known similar representations as double integral like
$$ \gamma~=~  \int \limits_0^1 \int \limits_0^1 \frac{1-x}{(1-x\, y ) (-\log x\,y)}\,dx\,dy.  $$
Is there a valid way to prove \eqref{1} by converting it into an integral representation ?
 A: $\newcommand\ga\gamma$As usual, let $H_n$ denote the $n$th harmonic number. Note that
$$\frac{s-k}{k (s n+k)}=\frac1n\frac1k-\frac{n+1}n\frac1{n s+k}$$
and hence
$$
\begin{aligned}
&\sum_{k,s=1}^n\frac{s-k}{k (s n+k)} \\
&=\frac1n\,\sum_{k,s=1}^n\frac1k-\frac{n+1}n\,\sum_{k,s=1}^n\frac1{n s+k} \\ 
&=H_n-\frac{n+1}n\,\sum_{s=1}^n\sum_{k=1}^n\frac1{n s+k} \\ 
&=H_n-\frac{n+1}n\,\sum_{s=1}^n (H_{n(s+1)}-H_{ns}) \\ 
&=H_n-\frac{n+1}n\,(H_{n(n+1)}-H_n) \\ 
&=\frac{2n+1}n\,H_n-\frac{n+1}n\,H_{n(n+1)}\\ 
&=\frac{2n+1}n\,(\ga+\ln n+o(1)) \\ 
&-\frac{n+1}n\,(\ga+\ln(n(n+1))+o(1)) \\ 
&=\ga+o(1)
\end{aligned}$$
(as $n\to\infty$), which proves your formula $(*)$.
Of course, $H_n$ can be represented as an integral; for instance, we have Euler's representation
$$H_n=\int _{0}^{1}{\frac {1-x^{n}}{1-x}}\,dx.$$

On the other hand, apparently, your Riemann-sum argument cannot work -- because the boundary effect here is nonnegligible. In particular, summing over the set $\{2,\dots,n\}\times\{2,\dots,n\}$ rather than over the set $\{1,\dots,n\}\times\{1,\dots,n\}$, we get
$$\lim_{n\to\infty}\sum_{k,s=\color{red}{\mathbf{2}}}^n\frac{s-k}{k (s n+k)}=\ga-1+\ln2\ne\ga.$$
