Riemann zeta for real argument between 0 and 1 using Mellin, with short asymptotic expansion

The Euler-Maclaurin formula gives the asymptotic expansion

\begin{equation} \sum_{k = 1}^{n} \frac{1}{k^{\sigma}} = \zeta(\sigma) + \sum_{k = 0}^{\infty} \frac{(-1)^k}{1-\sigma} \binom{1-\sigma}{k} B_k n^{1-\sigma-k}\,. \end{equation}

Thus $A = -B_1 = \frac{1}{2}$ and $B = -\frac{(-1)\sigma}{2}B_2 = \frac{\sigma}{12}$.

Taking a few more terms, we get

\begin{align} \zeta(\sigma) &= \sum_{k = 1}^n \frac{1}{k^{\sigma}} - \frac{n^{1-\sigma}}{1-\sigma} - \frac{1}{2n^{\sigma}} + \frac{\sigma}{12 n^{1+\sigma}} - \frac{\sigma(1+\sigma)(2+\sigma)}{720 n^{3+\sigma}} \\ &\qquad+ \frac{\sigma(1+\sigma)(2+\sigma)(3+\sigma)(4+\sigma)}{30240n^{5+\sigma}} + O(n^{-7-\sigma})\,. \end{align}

There is a standard technique that produces the complete asymptotic expansion for this sum and many others like it, which is to use harmonic sums and Mellin transforms.

Introduce the telescoping sum where $0\lt\sigma\lt 1$

$$S(x) = \sum_{k\ge 1} \left(\frac{1}{k^\sigma}- \frac{1}{(x+k)^\sigma}\right).$$

This sum has the property that

$$S(n) = \sum_{p=1}^n \frac{1}{p^\sigma},$$

so that $S(n)$ is the value we are looking for.

Re-write the sum as follows:

$$S(x) = \sum_{k\ge 1} \frac{1}{k^\sigma} \left(1-\frac{1}{(x/k+1)^\sigma}\right).$$

The sum term is harmonic and may be evaluated by inverting its Mellin transform.

Recall the harmonic sum identity $$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) = \left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$ where $g^*(s)$ is the Mellin transform of $g(x).$

In the present case we have $$\lambda_k = \frac{1}{k^\sigma}, \quad \mu_k = \frac{1}{k} \quad \text{and} \quad g(x) = 1 - \frac{1}{(1+x)^\sigma}.$$

It follows that $$\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1} \frac{1}{k^\sigma} \times k^s = \zeta(\sigma-s)$$ which has fundamental strip $\sigma-s > 1$ or $s < \sigma-1.$

We need the Mellin transform $g^*(s)$ of $g(x)$ which is $$\int_0^\infty \left(1 - \frac{1}{(1+x)^\sigma}\right) x^{s-1} dx$$

which is immediately seen to be a beta function integral with value $$g^*(s) = - \frac{1}{\Gamma(\sigma)} \Gamma(s)\Gamma(\sigma-s)$$

and fundamental strip $\langle -1, 0 \rangle.$ We check that the abscissa of convergence of the zeta function term is right in this strip as required. It follows that the Mellin transform $Q(s)$ of the harmonic sum $S(x)$ is given by

$$Q(s) = - \frac{1}{\Gamma(\sigma)} \Gamma(s)\Gamma(\sigma-s) \zeta(\sigma-s).$$

The Mellin inversion integral here is $$\frac{1}{2\pi i} \int_{\sigma-1-\varepsilon-i\infty} ^{\sigma-1-\varepsilon+i\infty} Q(s)/x^s ds$$

which we evaluate by shifting it to the right for an expansion at infinity.

First treat the pole from the zeta function term at $s=\sigma-1$, which has

$$\mathrm{Res}(Q(s)/x^s; s=\sigma-1) = -\frac{1}{\Gamma(\sigma)} \Gamma(\sigma-1)\Gamma(1) \times -1 \times x^{1-\sigma} = -\frac{1}{1-\sigma} x^{1-\sigma}.$$

For the pole at $s=0$ from the simple gamma function term we obtain

$$\mathrm{Res}(Q(s)/x^s; s=0) = -\frac{1}{\Gamma(\sigma)} \Gamma(\sigma) \zeta(\sigma) = -\zeta(\sigma).$$

For the pole at $s=\sigma$ from the compound gamma function term we obtain

$$\mathrm{Res}(Q(s)/x^s; s=\sigma) = -\frac{1}{\Gamma(\sigma)} \Gamma(\sigma) \times -1 \times\zeta(0) \times \frac{1}{x^\sigma} = -\frac{1}{2} \frac{1}{x^\sigma}.$$

The remaining poles are at $s = q+\sigma$ where $q\ge 1$ and contribute

$$\mathrm{Res}(Q(s)/x^s; s=q+\sigma) = -\frac{1}{\Gamma(\sigma)} \Gamma(\sigma+q) \frac{(-1)^{q+1}}{q!} \zeta(-q) \frac{1}{x^{q+\sigma}} \\ = - \prod_{p=0}^{q-1} (p+\sigma) \times \frac{(-1)^{q+1}}{q!} (-1)^q \frac{B_{q+1}}{q+1} \frac{1}{x^{q+\sigma}} = {q+\sigma-1\choose q} \frac{B_{q+1}}{q+1} \frac{1}{x^{q+\sigma}}.$$

The zero values of the Bernoulli numbers correctly represent cancelation of the gamma function poles by the trivial zeros of the zeta function.

Setting $x=n$ and observing that the shift to the right produces a minus sign we obtain the following asymptotic expansion:

$$S(n) = \sum_{p=1}^n \frac{1}{p^\sigma} \\ \sim \frac{1}{1-\sigma} n^{1-\sigma} + \zeta(\sigma) + \frac{1}{2} \frac{1}{n^\sigma} - \sum_{q\ge 1} {q+\sigma-1\choose q} \frac{B_{q+1}}{q+1} \frac{1}{n^{q+\sigma}}.$$

Actually computing the Bernoulli number terms we get for the example by OP with $\sigma=1/3$ the expansion

$$3/2\,{n}^{2/3}+\zeta \left( 1/3 \right) +1/2\,{\frac {1}{\sqrt [3]{n}}} -1/36\,{n}^{-4/3}+{\frac {7}{4860\,{n}^{10/3}}} \\-{\frac {13}{26244\,{n}^{16/3}}} +{\frac {247}{590490}{n}^{-{\frac {22}{3}}}} -{\frac {6175}{9565938}{n}^{-{\frac {28}{3}}}} \\+{\frac {406999}{258280326}{n}^{-{\frac {34}{3}}}} -{\frac {12966835}{2324522934}{n}^{-{\frac {40}{3}}}}+\cdots$$

This MSE link points to a series of similar calculations.

This is indeed true. In fact, this formula can be derived from the Euler summation formula with $f(x)=x^{-s}$ with $0<s<1$. We have

$$\begin{align*} \sum_{n\leq x}\frac{1}{x^s}&=\int_1^x\frac{dt}{t^s}-s\int\frac{t-[t]}{t^{s+1}}dt+1-\frac{x-[x]}{x^s} \\ &=\frac{x^{1-s}}{1-s}-\frac{1}{1-s}+1-s\int_1^\infty\frac{t-[t]}{t^{s+1}}dt+O(x^{-s}). \end{align*}$$

If we take the limit as $x\rightarrow\infty$ we find that

$$1-\frac{1}{1-s}-s\int_1^\infty\frac{t-[t]}{t^{s+1}}dt=\zeta(s)$$

so we then have

$$\lim_{x\rightarrow\infty}\left(\sum_{n\leq x}\frac{1}{n^s}-\frac{x^{1-s}}{1-s}\right)=\zeta(s),\;\;0<s<1$$

as desired.

Note that $[t]$ is the greatest integer function. This proof can be found in Introduction to Analytic Number Theory by Apostol.