Optimal quadrature rule for heavy tail measure

I'm looking for a well-thought quadrature rule for this measure $$d\mu(t)=\frac{dt}{t^s}$$ for $$s\in(0,1)$$, the underlying motivation is to compute this integral $$\lambda^{s-1}=\frac{1}{\Gamma(1-s)}\int_{0}^{\infty} e^{-t\lambda}\frac{dt}{t^s}$$ with $$\lambda > 0$$, see this reference, Page 5, equation (5), https://arxiv.org/abs/1808.05159. So I'm looking for a way of dealing with the singularity at $$t=0$$ as well as the limit to $$\infty$$. Observe that $$\int_{T}^\infty \frac{dt}{t^s} = \infty$$ I was using the QuadGK.jl package and it is working just fine

λ = 2; s = 0.8
λs = λ^(s-1)

f(t) = exp(-λ*t)/t^s/gamma(1-s)


The predicted error and the real error is the same order of magnitude. However, when using the package for another formula $$\lambda^s = \frac{1}{\Gamma(-s)}\int_0^\infty (e^{-t\lambda}-1)\frac{dt}{t^{1+s}}$$ same reference as before. I again do the quadrature

λs = λ^s
g(t) = (exp(-λ*t)-1)/t^(1+s)/gamma(-s)


it works just fine, but this time the real error is of order 1e-3 while the predicted error is around 1e-8. And if I put a positive matrix instead of a scalar

A = # some positive definite matrix
h(t) = (exp(-A.*t)-I)/t^(1+s)/gamma(-s)


I get an ERROR: DomainError with 0.5: integrand produced NaN in the interval (0.0, 1.0), and, to my understanding, quadgk does not evaluate the function at $0$ as per the documentation https://juliamath.github.io/QuadGK.jl/stable/. And if I write

res, err = quadgk(h,1e-50,Inf)


the integral blows up. If I augment the lower bound to 1e-20 then I get a norm2 error of 0.24 w.r.t the fractional matrix, if I keep augmenting the lower bound (1e-10, 1e-5), the norm2 error keeps augmenting. I computed the error using the exact eigenvalues and eigenvectors. Has somebody got an intuition to where the instability behind the scenes might be? I would also like to know whether the Gauss-Kronrod quadrature is optimal for this integral, and if I'm doing fine, as I'm relatively new to this :)

EDIT (following a comment): unfortunately does not make a lot of difference, I leave the fully reproducible code here

using QuadGK
using SpecialFunctions

s = 0.8
nx = 3; Δx = 1/(nx+1) # change nx to whatever
A = -[2. -1 0; -1 2 -1; 0 -1 2]/Δx^2
Χ = zeros(Float64,nx,nx)
for p=1:nx
Χ[:,p] = sin.(p*pi*(1:nx)*Δx)
end
Λ = diagm(2*(cos.((1:nx)*pi*Δx).-1)/Δx^2)
As = -Χ*(diagm(diag(-Λ).^s))/Χ
h(t) = (exp(A.*t)-I)/t^(1+s)/gamma(-s)
res2 = res21 + res22; err2 = err21 + err22
@show res2
@show norm(res2+As)/norm(As)
@show norm(res1+As)/norm(As)
@show norm(res1-res2)


You can split the integral up such that you don't get the singularity at $$t=0$$. Something like this

using QuadGK
using SpecialFunctions

λ = 2.0; s = 0.8
λs = λ^(s-1)

f(t) = exp(-λ*t)/t^s/gamma(1-s)

λs = λ^s
g(t) = (exp(-λ*t)-1)/t^(1+s)/gamma(-s)

A = [2 -1 0; -1 2 -1; 0 -1 2] # some positive definite matrix
h(t) = (exp(-A.*t)-I)/t^(1+s)/gamma(-s)

singularity_integral, singularity_error = quadgk(h, 1e-10, 1e-5)
regular_integral, regular_error = quadgk(h, 1e-5, Inf)
@show result = singularity_integral + regular_integral
@show error = singularity_error + regular_error


which will print

(res, err) = quadgk(f, 0, Inf) = (0.8705505365052336, 1.2164780075652479e-8)
(res, err) = quadgk(g, 0, Inf) = (1.7402976255082268, 2.515036776166618e-8)
result = singularity_integral + regular_integral = [1.6837899364901885 -0.7050512093960599 -0.03988519645256732; -0.7050512093960599 1.6439047399317526 -0.7050512098966073; -0.03988519645256732 -0.7050512098966073 1.683789936502625]
error = singularity_error + regular_error = 3.822238569033589e-8


Alternatively, you can also look into something like Cauchy principal values and this example

--Edit: update to new information--

I'm afraid I can't reproduce the same error. When I run your code (I've added using LinearAlgebra for diagm) I get

res2 = [15.354655632950173 -6.419790221101635 -0.3665284746983265; -6.419790221101635 14.988127158828647 -6.4197902215468154; -0.3665284746983265 -6.4197902215468154 15.354655632887033]
norm(res2 + As) / norm(As) = 0.018647374162762107
norm(res1 + As) / norm(As) = 0.01864737561229984
norm(res1 - res2) = 4.9872522079524005e-8
4.9872522079524005e-8