Adding this as an answer because it's kind of long for a comment, but it does not answer your question completely:
First, is the lower bound supposed to be one or zero? I was able to evaluate a closed form solution for generic p >1 for at least just the integral w.r.t. j so at least one of the integrations can be done:
ClearAll["Global`*"];
$Assumptions = {d > 0, s > 0, p > 1};
h[i_] = (i + 1)^-p;
f[i_, s_] = h[i] Exp[-2 s h[i]];
integrand = (f[i, s] - f[j, s])^2;
justJ = FullSimplify[Integrate[integrand, {j, 0, d}]]
(*(1/(16 p^3 s^2))(p (16 E^(-2 (1 + (1 + i)^-p) s) (1 + i)^-p p s -
16 (1 + d) E^(-2 ((1 + d)^-p + (1 + i)^-p) s) (1 + i)^-p p s +
16 d E^(-4 (1 + i)^-p s) (1 + i)^(-2 p) p^2 s^2 -
E^(-4 s) (-1 + p + 4 p s) + (1 + d) E^(-4 (1 + d)^-p s) (-1 +
p + 4 (1 + d)^-p p s)) -
16 E^(-2 (1 + i)^-p s) (1 + i)^-p p s ExpIntegralE[1 + 1/p,
2 s] + (-1 + p) ExpIntegralE[1 + 1/p, 4 s] +
16 (1 + d) E^(-2 (1 + i)^-p s) (1 + i)^-p p s ExpIntegralE[1 + 1/p,
2 (1 + d)^-p s] - (1 + d) (-1 + p) ExpIntegralE[1 + 1/p,
4 (1 + d)^-p s])*)
This took about 10s on my computer.
Edit: It takes only another ~15s to evaluate:
final = FullSimplify[Integrate[justJ, {i, 0, d}]]
(* (1/(p^2))(1 +
d)^(-2 p) (-2 d (1 + d)^(2 p) p ExpIntegralE[-1 + 1/p, 4 s] +
2 d (1 + d) p ExpIntegralE[-1 + 1/p, 4 (1 + d)^-p s] -
2 ((1 + d)^
p ExpIntegralE[1/p, 2 s] - (1 + d) ExpIntegralE[1/p,
2 (1 + d)^-p s])^2) *)
But I cannot get Limit[final/d, d -> Infinity]
to evaluate. I am also not 100% sure if I've made any mathematical sin by doing the integrals first and then taking the limit as d->Infinity
Final Edit:
after going to get coffee, (so < 5 min), the limit is evaluated as:
Limit[final/d, d -> Infinity]
(*(2^(-3 + 2/p) s^(-2 + 1/p) (Gamma[2 - 1/p] - Gamma[2 - 1/p, 4 s]))/p*)
sorry for all of the iterative edits on this response.