## Abstract

We approximate each generalized Stieltjes constant *γ*_{n}(*a*) by means of a finite sum involving Bernoulli numbers. Such an approximation has a quasi-geometric rate of convergence, which improves as Re(*a*) increases. A more detailed analysis, including numerical computations, is carried out for the constants *γ*_{0}(1) and *γ*_{1}(1). The key point in the proof is a probabilistic representation of the aforementioned constants, obtained as a consequence of a differential calculus concerning the gamma process.

## 1. Introduction and main result

Let
be the Hurwitz zeta function. The Laurent expansion of *ζ*(*s*,*a*) about its simple pole at *s*=1 is given by
1.1where the constants (*γ*_{n}(*a*))_{n≥0} are known as generalized Stieltjes constants. For *a*=1, *ζ*(*s*,1)=*ζ*(*s*) is the Riemann zeta function and the constants *γ*_{n}(1)=*γ*_{n} are called Stieltjes or generalized Euler constants, since *γ*_{0} is the Euler–Mascheroni constant. For any *n*=1,2,…, Berndt (1972) showed that
1.2thus extending a result pointed out by Stieltjes (1905) for *γ*_{n−1}. Different integral or series representations for *γ*_{n}(*a*) can be found in Israilov (1979, 1981), Zhang & Williams (1994), Coffey (2006, 2009*a*, 2011) and the references therein.

One of the main questions concerning generalized Stieltjes constants is to obtain upper bounds for |*γ*_{n}(*a*)| useful for large values of *n*. In this regard, we mention the explicit estimates provided by Zhang & Williams (1994) and Israilov (1979, 1981), the last ones in terms of Bernoulli numbers. As far as we know, the best explicit upper bound to date for |*γ*_{n}|, *n*≥10 is
1.3where and ⌊*x*⌋ stands for the integer part of *x*. Such a bound follows by gathering together the results by Matsuoka (1985) and Adell (2011). The order of magnitude of the right-hand side in (1.3) is
In two recent papers, Knessl & Coffey (2011*a*,*b*) have obtained the leading asymptotic form of *γ*_{n}(*a*). More precisely, denote by *f*(*n*)∼*g*(*n*) whenever *f*(*n*)/*g*(*n*)→1, as . After some computations, these two authors show that
1.4

A different, but closely related problem, is to estimate *γ*_{n}(*a*) for each *n*=0,1… and Re(*a*)>0. In this respect, the rate of convergence in (1.2) is rather slow (see the comments following theorem 1.1). From a numerical point of view, Kreminski (2003) has computed *γ*_{n}(*a*) for *n*=0,1,…,3200 (see also Choudhury 1995). Coffey (2007) has obtained rapidly convergent expressions for *γ*_{n} in terms of Bernoulli numbers. Also, series representations for *γ*_{0} and *γ*_{1} with an exponential rate of decay can be found in Coffey (2009*a*,*b*). Finally, Coffey (2011) has provided hypergeometric summation representations for *γ*_{1} and *γ*_{2} with summands of the order of *O*(*n*^{−6}) and *O*(*n*^{−3}), respectively.

The aim of this paper is to approximate *γ*_{n}(*a*), for each *n*=0,1,… and Re(*a*)>0, by means of a finite sum, with a quasi-geometric rate of convergence. Let us describe the ingredients appearing in such an approximating sum. In first place, we consider the function
1.5where (*B*_{k})_{k≥0} is the sequence of Bernoulli numbers, i.e.
1.6On the other hand, for any *n*=1,2,…, *k*=0,1,…, and Re(*s*)>−1, we consider the coefficients
1.7where the sum is extended to those non-negative integers *k*_{1},…,*k*_{n}, such that *k*_{1}+⋯+*k*_{n}=*k*, and it is understood that . A probabilistic meaning of such coefficients will be given in lemma 3.1 in §3. Observe, in particular, that
1.8

Denote by and by ⌈*x*⌉ the ceiling of *x*, i.e. the smallest integer not less than *x*. From now on, we fix *n*=1,2,… and with *σ*=Re(*a*)>0, denote by *α*=1.615395⋯ the solution to the equation
1.9and denote by *m*≥2 an integer such that
1.10The assumptions (1.9) and (1.10) satisfied by *α* and *m* will become clear in §4. We are in a position to state the main result of this paper.

### Theorem 1.1

*With the preceding notations, we have
*1.11*where the main term of the approximation τ*_{n,a}*(m) is given by
*1.12*and the upper bound b*_{n,a}*(m) is given by
*1.13

Concerning theorem 1.1, some remarks are in order. For fixed *n*=1,2,… and with Re(*a*)>0, it follows from (1.11) and (1.12) that
This was already noticed by Zhang & Williams (1994) and explains why the rate of convergence in (1.2) is rather poor, specially for big *n*.

On the other hand, the role played by the parameter *a* in the Hurwitz zeta function can be summarized as follows. Whenever |*Im*(*a*)| remains bounded, we can approximate *γ*_{n−1}(*a*) by the finite sum *τ*_{n,a}(*m*) with a quasi-geometric rate of convergence
1.14since e^{−α}<*α*/2*π*, as follows from (1.9). Such a rate of convergence improves as *σ*=Re(*a*) increases. The number of terms of the approximating sum *τ*_{n,a}(*m*) is 2*m*+⌈*σ*⌉, approximately (recall that *m*+⌈*σ*⌉≥*n*, thanks to assumption (1.10)). For large values of |Im(*a*)|, we find in (1.14) a constant of the order of |Im(*a*)|^{n+1}. This kind of dependence on |Im(*a*)| is not unusual. For instance, in many algorithms to compute the Riemann zeta function *ζ*(*s*), we also find in the corresponding upper bounds certain constants of the order of (cf. Borwein (2000) and Coffey (2009*b*) among others). Even if Im(*a*)=0, we have constants of the order of 2^{n+1} in the upper bound of *b*_{n,a}(*m*) given in (1.13). These facts are not surprising, in view of the growth of *γ*_{n}(*a*) as *n* increases (recall formula (1.4)).

Looking at the classical Stieltjes constants *γ*_{n}=*γ*_{n}(1), a possible interesting question is the following: how many terms *m*=*m*(*n*) in the finite sum *τ*_{n,1}(*m*) are needed to approximate *γ*_{n} in a reasonable way, uniformly as *n* tends to infinity? Observe that taking terms, theorem 1.1 gives us an approximation of *γ*_{n} of the order of

To illustrate theorem 1.1, it may be of interest (as suggested by one of the referees) to specify the main terms *τ*_{1,1}(*m*) and *τ*_{2,1}(*m*) approximating the constants *γ*_{0} and *γ*_{1}, respectively. This is done in the following.

### Corollary 1.2

*In the setting of theorem 1.1, we have*
1.15and
1.16

As follows from (1.13) and (1.14), the orders of magnitude of the upper bounds *b*_{1,1}(*m*) and *b*_{2,1}(*m*) are
respectively. However, table 1 shows that such upper bounds overestimate the real differences |*γ*_{0}−*τ*_{1,1}(*m*)| and |*γ*_{1}−*τ*_{2,1}(*m*)|, respectively (the computations given in table 1 have been performed with MAPLE v. 12.0). In the general context of theorem 1.1, we do not know if the rate of convergence in approximating *γ*_{n−1}(*a*) by the finite sum *τ*_{n,a}(*m*) is in fact geometric for each fixed *n*=1,2,… and with Re(*a*)>0.

The paper is organized as follows. In the following section, we give a probabilistic representation of *γ*_{n}(*a*), obtained as a consequence of the differential calculus for linear operators represented by stochastic processes, particularly gamma processes, developed in Adell & Lekuona (2000). Such a calculus has already found applications in analytical number theory (cf. Adell & Jodrá 2008) and information theory (cf. Adell *et al.* 2010), among other fields. Starting from the aforementioned probabilistic representation, it becomes clear how to distinguish between the main term of the approximation (evaluated in §3) and the remainder term (bounded in §4). Finally, the proof of corollary 1.2 is given in §5.

## 2. Probabilistic representation for *γ*_{n}(*a*)

From now on, the following notations and assumptions will be used. By , we denote the set of non-negative integers and by . Every function *ϕ* is a complex-valued function defined on or (occasionally, on ) for which the expectations under consideration are finite. The indicator function of the set *A* is denoted by 1_{A}. All of the random variables appearing in a same expression are supposed to be mutually independent.

Let (*X*_{t})_{t≥0} be a gamma process, i.e. a stochastic process starting at the origin, with independent stationary increments, right-continuous non-decreasing paths, and such that for each *t*>0, the random variable *X*_{t} has the gamma density
2.1Observe that
2.2

In order to state the differential calculus for the gamma process, let *U* and *T* be two random variables such that *U* is uniformly distributed on [0,1] and *T* has the exponential density *ρ*_{1}(*θ*). By (*U*_{k})_{k≥1} and (*T*_{k})_{k≥1}, we denote two sequences of independent copies of *U* and *T*, respectively. Set
2.3By Fubini’s theorem and (2.2), we have
2.4Thus, the Laplace transform of *X*_{t} can be written as
This means that the gamma process is a centred subordinator with characteristic random variable *T* (see Feller (1966, p. 450) or Adell & Lekuona (2000, §4.3.1)). It therefore follows from Adell & Lekuona (2000), (corollary 1 and proposition 4) that we have the Taylor expansion
2.5where *S*_{n} is defined in (2.3). In the following result, we give a probabilistic representation of *γ*_{n}(*a*) in terms of the function
2.6where the function *f*(*x*) is defined in (1.5).

### Proposition 2.1

*For any* *we have*
2.7

### Proof.

Let *t*>0. Using the gamma representation for the Hurwitz zeta function (see Ivić 1985, p. 99), we have
2.8On the other hand, the Taylor expansion in (2.5) for *ϕ*=*f*_{a} and *r*=0 gives us
since *f*_{a}(*X*_{0})=*f*_{a}(0)=1. The result follows by comparing this expansion with that in (1.1) for *s*−1=*t*, and taking into account (2.8). *

The basic idea to prove theorem 1.1 comes from representation (2.7). The direct application of this formula, however, poses several difficulties. The first one is that the derivatives of the function *f*_{a} become very involved for large *n*. The second is that such derivatives attain huge values near the origin. In fact, we have for even *n*
as follows from (1.5) and (1.6). Such large values cannot be compensated by the small probabilities , , for certain positive constants *c* and . Finally, the expansion in (1.5) is valid only for *x*∈[0,2*π*) and, therefore, a full expansion of *f*_{a} in terms of the Bernoulli numbers is not possible on .

For the previous reasons, we consider the finite sum
2.9and decompose the function *f*_{a}(*x*), for *x*≥0, into
2.10where
2.11and
2.12For notational simplicity, the dependence on *m* and *a* of such functions is dropped (recall the conventions before (1.9) and (1.10)). As we will see in proposition 4.5, we need *m*+⌈*σ*⌉−1 terms in the sum on the right-hand side in (2.9) in order to give an upper bound in theorem 1.1 having a homogeneous order of magnitude. The functions *g* and *h*, which are easy to derive, will give rise to the main term *τ*_{n,a}(*m*) in theorem 1.1, whereas the function *ψ*, which is smooth at infinity and at the origin, will account for the remainder term.

## 3. Evaluation of the main term

In view of representation (2.7) and the comments at the end of §2, the approximant *τ*_{n,a}(*m*) of *γ*_{n−1}(*a*) is obtained by computing the *n*th derivatives *g*^{(n)}(*x*) and *h*^{(n)}(*x*), replacing *x* by the random variable *S*_{n} defined in (2.3), and then taking expectations. In this last respect, the following auxiliary result will be needed.

### Lemma 3.1

*Let* *with* Re(*s*)>−1. *For any* *we have*
3.1

*In addition, for any* *we have*
3.2*where c*_{n}(*s*,*k*) *is defined in* (1.7).

### Proof.

Since *T* has the exponential density *ρ*_{1}(*θ*), it follows from (2.1) and (2.2) that
3.3Fix *T*=*t*>0. Applying repeatedly integration by parts, we obtain
3.4Recall that the random variables *U* and *T* are independent. Thus, replacing *t* by *T* in (3.4) and taking expectations, we have from (2.4) and (3.3)
which shows (3.1). Using (2.3), (3.1) and the independence of the random variables involved, we get
thus concluding the proof. ■

### Remark 3.2

Formulas (3.1) and (3.2) make sense for *s*=0. As follows from the proof of lemma 3.1, the values of the right-hand sides in (3.1) and (3.2) at *s*=0 are, respectively,
and

### Proposition 3.3

*Let g*(*x*) *be as in (2.11). For any* *we have*

### Proof.

Let *x*≥0. For any *k*=0,1,…,*m*−1, denote by
3.5Observe that
3.6

On the other hand, it follows from (1.8) and lemma 3.1 that 3.7as well as 3.8The result readily follows from (2.11) and (3.6)–(3.8). ■

## 4. Estimate of the remainder term

Having in mind proposition 2.1 and the decomposition of the function *f*_{a} given in (2.9)–(2.12), the remainder term of *γ*_{n−1}(*a*) is
We shall bound this quantity separately on the events {*S*_{n}≤*α*} and {*S*_{n}>*α*}, where *α* is defined in (1.9). Previously, an estimate for the tail series of the polygamma function with real argument will be needed. It is known that this function can be written as a finite sum involving the Stirling numbers of the second kind or the Eulerian numbers (see Gould (1978) and Wang *et al.* (2010)). However, the estimate in the following auxiliary result, obtained by easy arguments, will suffice for our purposes.

### Lemma 4.1

*Let q*∈(0,1) *and β*>0 *with q* e^{β}<1. *Let* *with j*/*β*≤*M*. *Then*,

### Proof.

The function *F*(*x*)=*x*^{j} e^{−βx}, *x*≥0 attains its maximum at *x*=*j*/*β* and decreases in the interval . Since *j*/*β*≤*M*, we have
which shows the result. ■

### Proof.

Let 0≤*x*≤*α*. From (2.9) and (2.12), we see that
4.1where *h*_{k}(*x*) is defined in (3.9). Since *k*≥*m*+⌈*σ*⌉≥*n* (recall assumption (1.10)), we have from (3.10) and (4.1)
4.2On the other hand, it follows from (1.6) that
4.3where the last inequality follows by setting *q*=*α*/2*π* and *β*=1 in lemma 4.1 (observe that the assumptions in this lemma are fulfilled, thanks to hypotheses (1.9) and (1.10)). Again by (1.9), we have
4.4We therefore conclude from (4.2)–(4.4) that
4.5Finally, we have as in (3.7)
4.6The conclusion follows from (4.5) and (4.6). ■

### Remark 4.3

The estimate in (4.6) could be improved by using Hölder’s inequality and a subsequent estimate of *P*(*S*_{n}≤*α*). However, the effort is not worthy, because the order of magnitude of the upper bound in proposition 4.2 is better than that in proposition 4.5 below.

Denote by 4.7

### Proof.

Let *x*>*α*. Observe that
4.8where *g*_{m+k}(*x*) is defined in (3.5). It follows from (3.6) that
4.9On the other hand, the function *r*(*x*)=*x* e^{−(m+k+σ−1)x}, *x*≥0, decreases for *x*≥1/(*m*+*k*+*σ*−1). By (1.9) and (1.10), we have *m*≥2, *σ*>0, and *α*=1.615395…. Thus, *m*+*k*+*σ*−1>1>*α*^{−1}. We therefore have from (4.9)
4.10Using (1.10) and the inequality 1+*x*≤e^{x}, *x*≥0, we have
4.11Thus, the conclusion follows from (4.8), (4.10) and (4.11). ■

In the following result, we shall use the well-known Poisson-gamma relation (cf. Johnson *et al.* 1993, p. 164), that is,
4.12

### Proof.

Let *x*>*α*. Recall that
4.13where *h*_{k}(*x*) is defined in (3.9). It follows from (3.10), (4.12) and the fact that *x*>*α*
This, together with (4.13) and (1.6), implies that
4.14The conclusion follows from (4.14) and (1.9). ■

### Proof of Theorem 1.1.

By proposition 2.1 and the decomposition of the function *f*_{a} given in (2.9)–(2.12), we have
4.15The first two expectations on the right-hand side in (4.15) have been computed in propositions 3.3 and 3.4, respectively, and constitute the main term of the approximation *τ*_{n,a}(*m*) defined in (1.12). The third expectation has been estimated in propositions 4.2, 4.4 and 4.5, and gives rise to the upper bound in *b*_{n,a}(*m*) defined in (1.13). The proof is complete. ■

## 5. Proof of corollary 1.2

Corollary 1.2 follows from (1.12), proposition 5.1 below and some simple computations.

### Proposition 5.1

*For any m*=1,2,… *and k*=2,3,…, *we have*
5.1*and*
5.2

## Acknowledgements

This work has been partially supported by research grants MTM2008-06281-C02-01/MTM, DGA E–64, UJA2009/12/07 (Universidad de Jaén and Caja Rural de Jaén), and by FEDER funds. The author thanks the referees for their helpful suggestions and comments, which greatly improved the final outcome.

- Received September 12, 2011.
- Accepted January 3, 2012.

- This journal is © 2012 The Royal Society