A Conjectured Integer Sequence Arising From the Exponential Integral

Richard P. Brent Australian National University

Canberra ACT 2601 Australia

ei@rpbrent.com

M. L. Glasser Clarkson University

Potsdam, NY 13699-5820

USA

lglasser@clarkson.edu

Anthony J. Guttmann School of Mathematics and Statistics

The University of Melbourne Vic. 3010 Australia

guttmann@unimelb.edu.au

Abstract

Let f_{0}(z) = exp(z/(1 − z)), f_{1}(z) = exp(1/(1 − z))E_{1}(1/(1 − z)), where E_{1}(x) = ∫_{x}^{∞} e^{−t}t^{−1} dt. Let a_{n} = [z^{n}]f_{0}(z) and b_{n} = [z^{n}]f_{1}(z) be the corresponding Maclaurin series coefficients. We show that a_{n} and b_{n} may be expressed in terms of confluent hypergeometric functions.

We consider the asymptotic behaviour of the sequences (a_{n}) and (b_{n}) as n → ∞, showing that they are closely related, and proving a conjecture of Bruno Salvy regarding (b_{n}).

Let ρ_{n} = a_{n}b_{n}, so ∑ ρ_{n}z^{n} = (f_{0} ⊙f_{1})(z) is a Hadamard product. We obtain an asymptotic expansion 2n^{3/2}ρ_{n} ∼ − ∑ d_{k}n^{−k} as n → ∞, where d_{k} ∈ ℚ, d_{0} = 1. We conjecture that 2^{6k}d_{k} ∈ ℤ. This has been verified for k ≤ 1000.

1 Introduction

We consider two analytic functions,

f_{0}(z) := e^{z/(1−z)} = e^{−1} e^{1/(1−z)}

and

f_{1}(z) := e^{x}E_{1}(x), where x := 1/(1 − z) and E_{1}(x) := ∫_{x}^{∞} \frac{e^{−t}}{t} dt.

These functions are regular in the open disk D = {z ∈ ℂ : |z| < 1}. We write their Maclaurin coefficients as a_{n} := [z^{n}]f_{0}(z) and b_{n} = [z^{n}]f_{1}(z). Thus, in the disk D, f_{0}(z) = ∑_{n≥0} a_{n}z^{n} and f_{1}(z) = ∑_{n≥0} b_{n}z^{n}.

The functions f_{0}(z) and f_{1}(z) satisfy the same third-order linear differential equation with polynomial coefficients. Thus, the sequences (a_{n}) and (b_{n}) are D-finite and satisfy the same recurrence relation (for sufficiently large n).

There are several entries in the On-Line Encyclopedia of Integer Sequences (OEIS) related to the rational sequence (a_{n})_{n≥0}. The numerators are OEIS \underline{A067764}, and the denominators are OEIS \underline{A067653}. The integers n!a_{n} are given by OEIS \underline{A000262} and, with alternating signs, by OEIS \underline{A293125}. The numbers (b_{n})_{n≥0} are unlikely to be rational.^{1}

The numbers a_{n} and b_{n} may be expressed in terms of confluent hypergeometric functions. If M(a, b, z) = _{1}F_{1}(a; b; z) and U(a, b, z) are standard solutions of Kummer’s differential equation, then Lemmas 1–2 show that a_{n} = e^{−1}M(n + 1, 2, 1) and b_{n} = −Γ(n)U(n, 0, 1).

We are interested in the asymptotics of a_{n} and b_{n} for large n. Perron [17], following Fejér [8], showed that

a_{n} ∼ \frac{e^{2\sqrt{n}}}{2n^{3/4}\sqrt{πe}} .

Salvy^{2} conjectured that b_{n} is of order e^{−2\sqrt{n}}n^{−3/4}. We have verified this conjecture. In fact,

b_{n} ∼ −\frac{\sqrt{πe}}{n^{3/4}e^{2\sqrt{n}}} .

1In particular, b_{0} = G, where G := eE_{1}(1) ≈ 0.596 is the Euler-Gompertz constant, whose decimal digits are given by OEIS \underline{A073003}. We have b_{n} = a_{n}G − a′_{n}, where a′_{n} ∈ ℚ and a′_{n} satisfies essentially the same recurrence as a_{n}, but with different initial conditions. Clearly b_{n} ∈ ℚ if and only if G ∈ ℚ. All that is known is that at least one of γ and G is irrational [1, 18].

2Bruno Salvy, email to A. J. Guttmann et al., May 28, 2018.

A function of the form f(n) = exp(αn^{θ+o(1)}) for α ≠ 0, θ ∈ (0, 1), is called a stretched exponential in the physics/statistics literature (the term sub-exponential is used in complexity theory). Thus, a_{n} and b_{n} are stretched exponentials, with α = ±2 and θ = 1/2.

The motivation for this paper stems from some enumeration problems in algebraic combinatorics and mathematical physics. Many such problems involve ordinary generating functions of power series A(x) = ∑_{n≥0} A_{n}x^{n} in which A_{n} ∼ cμ^{n}n^{g}. In such cases, assuming that g is not a negative integer, one can write

A(x) ∼ c Γ(1 + g)(1 − μx)^{−(1+g)}

as x → 1/μ. However, in recent years there have been a number of examples, such as Av(1324) pattern-avoiding permutations [5], interacting partially-directed self-avoiding walks [16], and Dyck paths enumerated by maximum height [12], in which the corresponding generating function has coefficients behaving as B_{n} ∼ cμ^{n} exp(αn^{θ})n^{g}, with α < 0. The question then arises as to the asymptotic form of the generating function. The coefficients b_{n} considered in this paper are of the form just described, with θ = 1/2, and the underlying generating function is found. Corresponding results for other values of θ remain to be discussed.

Theorem 5 gives complete asymptotic expansions of a_{n} and b_{n}. These may be written as

a_{n} = \frac{F (n^{1/2})}{2n^{3/4}\sqrt{πe}} and b_{n} = −\frac{\sqrt{πe}}{n^{3/4}} F (−n^{1/2}),

where F (x) ∼ e^{2x} ∑_{k≥0} c_{k}x^{−k}, for certain constants c_{k} ∈ ℚ, c_{0} = 1. The c_{k} may be computed using Theorem 5 or Lemma 7.

The Hadamard product f_{0} ⊙f_{1} of f_{0} and f_{1} is the analytic function defined for z ∈ D by

(f_{0} ⊙f_{1})(z) = ∑_{n≥0} a_{n}b_{n}z^{n}.

The asymptotic expansions of a_{n} and b_{n} imply an asymptotic expansion for ρ_{n} := a_{n}b_{n} of the form

ρ_{n} ∼ −\frac{1}{2n^{3/2}} ∑_{k≥0} d_{k}n^{−k},

where d_{k} ∈ ℚ, d_{0} = 1 (see Corollary 9).

A dyadic rational is a rational number of the form p/q, where q is a power of two. Let Q_{2} := {j/2^{k} : j, k ∈ ℤ} denote the set of dyadic rationals.

We conjecture, from numerical evidence for k ≤ 1000, that d_{k} ∈ Q_{2}. More precisely, defining r_{k} := 2^{6k}d_{k}, Conjecture 10 is that r_{k} ∈ ℤ. Remark 11 gives numerical evidence for a slightly stronger conjecture. In Theorem 17 we prove the weaker (but still nontrivial) result that k!r_{k} ∈ ℤ.

In Remark 13 we mention an analogous (easily proved) result for modified Bessel functions, where the product I_{ν}(x)K_{ν}(x) for fixed ν ∈ ℤ has an asymptotic expansion whose coefficients are in ℚ_{2}.

The connection with confluent hypergeometric (Kummer) functions is discussed in §2, and asymptotic expansions for a_{n} and b_{n} are considered in §3. In §4 we mention various recurrence relations, continued fractions, and closed-form expressions related to a_{n} and b_{n}. Finally, in §§5–6, we consider Hadamard products and discuss the conjecture mentioned above.

Some comments on notation: f(x) ∼ ∑_{k≥0} f_{k}x^{−k} means that the sum on the right is an asymptotic series for f(x) in the sense of Poincaré. Thus, for any fixed m > 0, f(x) = ∑_{k=0}^{m−1} f_{k}x^{−k} + O(x^{−m}) as x → ∞. The letters j, k, m, n always denote integers (except for n in Remark 4). The notation (x)_{n} for n ≥ 0 denotes the ascending factorial or Pochhammer symbol, defined by (x)_{n} := x(x + 1) ⋯ (x + n − 1).

2 Connection with hypergeometric functions

The numbers a_{n} and b_{n} may be expressed in terms of confluent hypergeometric functions (Kummer functions), for which we refer to [15, §13.2]. If M(a, b, z) and U(a, b, z) are standard solutions w(a, b, z) of Kummer’s differential equation zw′′ + (b − z)w′ − aw = 0, then Lemmas 1–2 below express a_{n} and b_{n} in terms of M(n + 1, 2, 1) and U(n, 0, 1).

Kummer [13] considered

M(a, b, z) = _{1}F_{1}(a; b; z) = ∑_{k≥0} \frac{(a)_{k} z^{k}}{(b)_{k} k!} ,  (1)

which is undefined if b is zero or a negative integer. In the case a ≠ b = 0, we can use the solution

zM(a + 1, 2, z) = lim_{b→0} \frac{b}{a}M(a, b, z).

Tricomi [22] introduced the function U(a, b, z) as a second (minimal) solution of Kummer’s differential equation. For our purposes it is convenient to use the integral representation [15, (13.4.4)] (valid for ℜ(a) > 0, ℜ(z) > 0)

U(a, b, z) = \frac{1}{Γ(a)} ∫_{0}^{∞} e^{−zt} t^{a−1} (1 + t)^{b−a−1} dt.  (2)

We remark that the functions M and U satisfy recurrence relations, known as “connection formulas”. For example, we mention [15, (13.3.1) and (13.3.7)], both (essentially) due to Gauss (see Erdélyi [7, §6.4 and §6.6]):

(b − a)M(a − 1, b, z) + (2a − b + z)M(a, b, z) − aM(a + 1, b, z) = 0,  (3)
U(a − 1, b, z) + (b − 2a − z)U(a, b, z) + a(a − b + 1)U(a + 1, b, z) = 0.  (4)

Lemmas 1–2 express a_{n} and b_{n} in terms of the Kummer functions M and U, respectively. Lemma 1 was stated, without proof, by Covo [6].

Lemma 1. If n ∈ ℤ, n ≥ 1, and a_{n} is as above, then

a_{n} = e^{−1}M(n + 1, 2, 1).  (5)

Proof. If we put a = n + 1, b = 2, and z = 1 in the connection formula (3), we see that ã_{n} := e^{−1}M(n + 1, 2, 1) satisfies the same recurrence (21) as a_{n}. Thus, to show that a_{n} = ã_{n} for all n ≥ 1, it is sufficient to show that a_{n} = ã_{n} for n ∈ {1, 2}. Now

ã_{1} = e^{−1}M(2, 2, 1) = e^{−1} ∑_{k≥0} \frac{(2)_{k}}{(2)_{k} k!} = 1 = a_{1},

and, similarly,

ã_{2} = e^{−1}M(3, 2, 1) = e^{−1} ∑_{k≥0} \frac{(3)_{k}}{(2)_{k} k!} = e^{−1} ∑_{k≥0} \frac{k + 2}{2 k!} = 3/2 = a_{2},

so the result follows.

Lemma 2. If n ∈ ℤ, n ≥ 1, and b_{n} is as above, then

b_{n} = −Γ(n) U(n, 0, 1).  (6)

Proof. We start with [15, (6.7.1)]:

I(a, b) := ∫_{0}^{∞} \frac{e^{−at}}{t + b} dt = e^{ab}E_{1}(ab), a, b > 0.

Note that, by definition, b_{n} = [z^{n}]I(1, 1/(1 − z)). Setting a = 1, b = 1/(1 − z), the term 1/(t + b) inside the integral can be rearranged as follows:

(t + \frac{1}{1 − z} )^{−1} = \frac{1 − z}{1 + t − tz} = \frac{1}{1 + t} − \frac{1}{t(1 + t)} (\frac{1}{1 − zt/(1 + t)} − 1) ,

and making the substitution s = t/(1 + t) gives

I(1, 1/(1 − z)) = ∫_{0}^{∞} \frac{e^{−t}}{1 + t} dt − ∫_{0}^{1} e^{−s/(1−s)} (\frac{z}{1 − zs}) ds = ∑_{n≥0} b_{n} z^{n}.

Thus, b_{0} = eE_{1}(1) and, for n > 0,

b_{n} = − ∫_{0}^{1} e^{−s/(1−s)} s^{n−1} ds.  (7)

Writing e^{−s/(1−s)} = e^{1−1/(1−s)} gives, for n > 0,

b_{n} = −e ∫_{0}^{1} e^{−1/(1−s)} s^{n−1} ds.  (8)

Substitute t = s/(1 − s) in (2), giving

Γ(a)U(a, b, z) = e^{z} ∫_{0}^{1} e^{−z/(1−s)} s^{a−1} (1 − s)^{−b} ds.  (9)

Comparison of (8) and (9) now gives b_{n} = −Γ(n) U(n, 0, 1).

Remark 3. We could prove Lemma 2 in the same manner as Lemma 1, using the connection formula (4) instead of (3), and the recurrence (23) instead of (21), but in order to verify the initial conditions we would have to resort to some explicit representation for U, such as the integral representation (2), so the proof would be no simpler.

Remark 4. We can generalize our definitions of a_{n} and b_{n} to permit n ∈ ℂ, using Lemmas 1– 2. Such generalizations do not seem particularly useful, so in what follows we continue to assume that n ∈ ℤ.

3 Asymptotic expansions of a_{n} and b_{n}

Theorem 5 gives the complete asymptotic expansions of a_{n} and b_{n} in ascending powers of n^{−1/2}. Wright [25] proved the existence of an asymptotic expansion of the form (10) for a_{n}, but did not state an explicit formula or algorithm for computing the constants c_{m} occurring in the expansion. For a more “algorithmic” approach, see Wyman [26].

Theorem 5. For positive integer n, if a_{n} and b_{n} are as above, then

a_{n} ∼ \frac{e^{2\sqrt{n}}}{2n^{3/4}\sqrt{πe}} ∑_{m≥0} c_{m}n^{−m/2}  (10)

and

b_{n} ∼ − \frac{\sqrt{πe}}{n^{3/4}e^{2\sqrt{n}}} ∑_{m≥0}(−1)^{m}c_{m}n^{−m/2},  (11)

where

c_{m} = (−1)^{m} ∑_{j=0}^{m} [h^{m−j}] exp(μ(h)) \frac{(m − 2j + 3/2)_{2j}}{4^{j}j!}  (12)

and

μ(h) = h^{−1} − (e^{h} − 1)^{−1} − \frac{1}{2} .  (13)

Remark 6. The function μ(h) defined by (13) could also be defined using Bernoulli numbers, since

μ(h) = − ∑_{k=1}^{∞} \frac{B_{2k}}{(2k)!}h^{2k−1} = −\frac{h}{12} + \frac{h^{3}}{720} − O(h^{5}).  (14)

The function exp(μ(h)) occurring in (12) has the Maclaurin expansion

exp(μ(h)) = 1 − \frac{h}{12} + \frac{h^{2}}{288} + \frac{67h^{3}}{51840} + O(h^{4}).  (15)

The numerators and denominators of the coefficients [h^{n}] exp(μ(h)) have been added to the OEIS as \underline{A321937} and \underline{A321938}, respectively.

Proof of Thm. 5. We first prove (11). From Lemma 2, b_{n} = −Γ(n) U(n, 0, 1). Temme [21, Sec. 3] gives a general asymptotic result for U(a, b, z^{2}) as a → ∞. We state Temme’s result for the case (a, b, z) = (n, 0, 1), which is what we need. Let c′_{k} := [h^{k}] exp(μ(h)). (Temme uses c_{k}, but this conflicts with our notation.) From Temme [21, (3.8)–(3.10)], we have

U(n, 0, 1) ∼ \frac{\sqrt{e}}{Γ(n)} ∑_{k≥0} c′_{k}Φ_{k}(n),  (16)

where

Φ_{k}(n) = 2n^{−(k+1)/2}K_{k+1}(2n^{1/2}),

and K_{ν} denotes the usual modified Bessel function.

From [15, (10.40.2)], K_{ν}(z) has an asymptotic expansion

K_{ν}(z) ∼ e^{−z}\sqrt{\frac{π}{2z}} ∑_{j≥0} \frac{(ν − j + 1/2)_{2j}}{j! (2z)^{j}} .  (17)

Setting ν = k and z = 2n^{1/2} in (17), we obtain

Φ_{k−1}(n) = 2n^{−k/2}K_{k}(2n^{1/2}) ∼ \frac{\sqrt{π}e^{−2\sqrt{n}}}{n^{1/4}} ∑_{j≥0} \frac{(k − j + 1/2)_{2j}}{j! 4^{j} n^{(j+k)/2}} .

Substituting this expression into (16), and grouping like powers of n, we obtain

b_{n} = −Γ(n) U(n, 0, 1) ∼ −\frac{\sqrt{πe}}{n^{3/4}e^{2\sqrt{n}}} ∑_{m≥0} ∑_{j=0}^{m} \frac{c′_{m−j} (m − 2j + 3/2)_{2j}}{j! 4^{j} n^{m/2}} .

Now, comparison with (11) shows that

(−1)^{m}c_{m} = ∑_{j=0}^{m} \frac{c′_{m−j} (m − 2j + 3/2)_{2j}}{j! 4^{j}} ,

which completes the proof of (11).

The proof of (10) is similar. We use Lemma 1 instead of Lemma 2, and Temme’s asymptotic result [21, (3.29)] for M(a, b, z^{2}) as a → ∞ instead of (16); the modified Bessel function I_{ν} replaces K_{ν}. From [15, (10.40.1)], I_{ν}(z) has an asymptotic expansion

I_{ν}(z) ∼ \frac{e^{z}}{\sqrt{2πz}} ∑_{j≥0}(−1)^{j} \frac{(ν − j + 1/2)_{2j}}{j! (2z)^{j}} ,  (18)

which replaces (17).

Theorem 5 gives an expression for c_{m} which (indirectly) involves Bernoulli numbers, in view of (14). Lemma 7 gives a different expression for c_{m} that is recursive, as the expression for c_{m} depends on the values of c_{j} for j < m, but has the advantage of avoiding reference to Bernoulli numbers. The idea of the proof is similar to that used in the “method of Frobenius” [10]. We omit the details, which may be found in [3, pp. 10–11].

Lemma 7. We have c_{0} = 1 and, for all m ≥ 1,

mc_{m} = [h^{m+3}] ∑_{j=0}^{m−1} c_{j}h^{j} ∑_{s∈\{±1\}}(1 + sh^{2})^{\frac{1−2j}{4}} exp (\frac{2}{h} ((1 + sh^{2})^{\frac{1}{2}} − 1)) .  (19)

Remark 8. Computation using (12) and, as a check, (19), gives

(c_{k})_{k≥0} = (1, −\frac{5}{48}, −\frac{479}{4608}, −\frac{15313}{3317760}, \frac{710401}{127401984}, −\frac{3532731539}{214035333120}, …) .

The numerators and denominators have been added to the OEIS as \underline{A321939} and \underline{A321940}, respectively. With the exception of c_{0} and c_{4}, the c_{k} all appear to be negative. This has been verified numerically for k ≤ 1000.

4 The Maclaurin coefficients a_{n} and b_{n}

The function f_{0}(z) is the exponential generating function counting several combinatorial objects, such as the number of “sets of lists”, i.e., the number of partitions of {1, 2, …, n} into ordered subsets, see Wallner [24, §5.3].

Observe that f_{0}(z) satisfies the differential equation

(1 − z)^{2}f′_{0}(z) − f_{0}(z) = 0,  (20)

and from this it is easy to see that the a_{n} satisfy a three-term recurrence

na_{n} − (2n − 1)a_{n−1} + (n − 2)a_{n−2} = 0 for n ≥ 2.  (21)

The initial conditions are a_{0} = a_{1} = 1. Thus

(a_{n})_{n≥0} = (1, 1, 3/2, 13/6, 73/24, 167/40, …).

The recurrence (21) holds for n ≥ 0 provided that we define a_{n} = 0 for n < 0. A closed-form expression, valid for n ≥ 1 (but not for n = 0), is

a_{n} = ∑_{k=1}^{n} \frac{1}{k!} \binom{n − 1}{k − 1}.

The constants a_{n} may be expressed in terms of the generalized Laguerre polynomials L_{n}^{(α)}(x) which, from [15, (18.12.13)], have a generating function

∑_{n≥0} z^{n}L_{n}^{(α)}(x) = (1 − z)^{−(α+1)}e^{−xz/(1−z)}.

With α = x = −1 we obtain ∑_{n≥0} z^{n}L_{n}^{(−1)}(−1) = e^{z/(1−z)}, so a_{n} = L_{n}^{(−1)}(−1).

Using the chain rule and the definition of f_{1}(z) in §1, we see that f_{1}(z) satisfies the differential equation

(1 − z)^{2}f′_{1}(z) − f_{1}(z) = z − 1,  (22)

which differs from (20) only in the right-hand side z − 1. Differentiating twice more with respect to z, we see that f_{0}(z) and f_{1}(z) both satisfy the same third-order differential equation

(1 − z)^{2}f′′′ + (4z − 5)f′′ + 2f′ = 0.

From (22), the b_{n} satisfy a recurrence

nb_{n} − (2n − 1)b_{n−1} + (n − 2)b_{n−2} = \begin{cases}1, & if n = 2; \\ 0, & if n ≥ 3. \end{cases}  (23)

This is essentially (i.e., for n ≥ 3) the same recurrence as (21), but the initial conditions b_{0} = G, b_{1} = G − 1 are different. Here G := eE_{1}(1) ≈ 0.596 is the Euler-Gompertz constant [14, §2.5].

We remark that computation of the b_{n} using the recurrence (23) in the forward direction is numerically unstable. A stable method of computation is to use an adaptation of Miller’s algorithm, originally used to compute Bessel functions. See Gautschi [11, §3] and Temme [20, §4].

As noted in §1, the b_{n} may be expressed as a_{n}G−a′_{n}, where a_{n} is as above, and a′_{n} satisfies essentially the same recurrence with different initial conditions. In fact,

na′_{n} − (2n − 1)a′_{n−1} + (n − 2)a′_{n−2} = \begin{cases}−1, & if n = 2; \\ 0, & if n ≥ 3. \end{cases}

The initial conditions are a′_{0} = 0, a′_{1} = 1. Thus

(a′_{n})_{n≥0} = (0, 1, 1, 4/3, 11/6, 5/2, 121/36, …).

From (11), b_{n} → 0 as n → ∞, so the sequence (a′_{n}/a_{n})_{n≥1} is a convergent sequence of rational approximations to G. The sequence of approximants is (1, 2/3, 8/13, 44/73, 100/167, …).

Bala [2] gives the continued fraction

1 − G = 1/(3 − 2/(5 − 6/(7 − ⋯ − n(n + 1)/(2n + 3) − ⋯))),

with convergents 1/3, 5/13, 28/73, 201/501, etc. The corresponding convergents to G are 2/3, 8/13, 45/73, 100/167, etc. We see that the n-th convergent is just a′_{n+1}/a_{n+1}. Theorem 5 implies that

G − a′_{n}/a_{n} = b_{n}/a_{n} ∼ −2πe^{1−4\sqrt{n}} as n → ∞.

We have contributed the sequence (n!a′_{n})_{n≥1} to the OEIS as \underline{A321942}.

5 The Hadamard product of f_{0} and f_{1}

Define ρ_{n} := a_{n}b_{n}. Thus ∑_{n=0}^{∞} ρ_{n}z^{n} is the Hadamard product (f_{0} ⊙f_{1})(z). From Lemmas 1– 2, we have

ρ_{n} = −e^{−1}Γ(n)M(n + 1, 2, 1)U(n, 0, 1).

Using Theorem 5, we can obtain a complete asymptotic expansion for ρ_{n} in decreasing powers of n. This is given in Corollary 9.

Corollary 9. We have

ρ_{n} ∼ − \frac{1}{2n^{3/2}} ∑_{k≥0} d_{k}n^{−k},

where

d_{k} = ∑_{j=0}^{2k}(−1)^{j}c_{j}c_{2k−j},

and c_{0}, …, c_{2k} are as in Theorem 5.

A computation shows that

(d_{k})_{k≥0} = (1, −7/32, 43/2048, −915/65536, …).

We observe that the d_{k} appear to be dyadic rationals More precisely, it appears that 2^{6k}d_{k} ∈ ℤ. Define a scaled sequence (r_{k})_{k≥0} by r_{k} := 2^{6k}d_{k}. Computation gives

(r_{k})_{k≥0} = (1, −14, 86, −3660, −1042202, −247948260, −108448540420, …).

This leads naturally to the following conjecture.

Conjecture 10. For all k ≥ 0, r_{k} ∈ ℤ.

The sequence of numerators of r_{k} has been added to the OEIS as \underline{A321941}. If Conjecture 10 holds, then the denominators are all 1, i.e., the denominators are given by \underline{A000012}. Remark 11. Conjecture 10 has been verified for all k ≤ 1000. We also showed numerically, for 3 ≤ k ≤ 1000, that r_{k} < 0 and r_{k} ≡ \binom{2k}{k} (mod 32). Remark 12. A problem that is superficially similar to our conjecture was solved by Tulyakov [23]. However, we do not see how to adapt his method to prove our conjecture. Remark 13. Corollary 9 is reminiscent of the result

I_{0}(x)K_{0}(x) ∼ \frac{1}{2x} ∑_{k≥0} e_{k,0} x^{−2k}

in the theory of Bessel functions [4, (1.2)]. The coefficients e_{k,0} are given by

e_{k,0} = \frac{(2k)!^{3}}{2^{6k}k!^{4}} ,

so 2^{4k}e_{k,0} ∈ ℤ. The modified Bessel functions I_{0}(x) and K_{0}(x) are solutions of the same ordinary differential equation xy′′ + y′ − xy = 0, but I_{0}(x) increases with x while K_{0}(x) decreases. This is analogous to the behaviour of a_{n}, which increases as n → ∞, and |b_{n}|, which decreases as n → ∞.

More generally, from [15, (10.40.6)], we have

I_{ν}(x)K_{ν}(x) ∼ \frac{1}{2x} ∑_{k≥0} e_{k,ν}x^{−2k},

where

e_{k,ν} = (−1)^{k}2^{−2k}(ν − k + 1/2)_{2k}\binom{2k}{k},

and 2^{4k}e_{k,ν} ∈ ℤ for ν ∈ ℤ.

6 Other expressions for d_{n}

Since (a_{n}) and (b_{n}) are D-finite, it follows that (ρ_{n}) is D-finite.^{3} In fact, ρ_{n} satisfies the 4-term recurrence

n^{2}(n − 1)(2n − 3)ρ_{n} = (n − 1)(2n − 1)(3n^{2} − 5n + 1)ρ_{n−1}
− (n − 2)(2n − 3)(3n^{2} − 5n + 1)ρ_{n−2}
+ (n − 2)(n − 3)^{2}(2n − 1)ρ_{n−3}  (24)

for n ≥ 3, with initial conditions ρ_{0} = G, ρ_{1} = G − 1, ρ_{2} = (9G − 6)/4.

The recurrence (24) can be simplified by defining σ_{n} := nρ_{n}. Then σ_{n} satisfies the slightly simpler recurrence

n(n − 1)(2n − 3)σ_{n} = (2n − 1)(3n^{2} − 5n + 1)σ_{n−1}
− (2n − 3)(3n^{2} − 5n + 1)σ_{n−2} + (n − 2)(n − 3)(2n − 1)σ_{n−3}  (25)

for n ≥ 3, with initial conditions σ_{0} = 0, σ_{1} = G − 1, σ_{2} = 9G/2 − 3. Also, Corollary 9 gives an asymptotic series for σ_{n}:

σ_{n} ∼ − \frac{1}{2n^{1/2}} ∑_{k≥0} d_{k}n^{−k}.  (26)

Using (25), we can give a recursive algorithm for computing the sequence (d_{n}) (and hence (r_{n})) directly, without computing the sequence (c_{n}).

3See Flajolet and Sedgewick [9, Appendix B.4], and Stanley [19, Theorem 2.10], for relevant background on D-finite sequences.

Lemma 14. We have d_{0} = 1 and, for all k ≥ 1,

8kd_{k} = − [h^{k+2}] ( ∑_{j=0}^{k−1} d_{j}h^{j}(B(h)(1 − h)^{−(j+1/2)}

+ C(h)(1 − 2h)^{−(j+1/2)} + D(h)(1 − 3h)^{−(j+1/2)})),  (27)

where

B(h) = −6 + 13h − 7h^{2} + h^{3} = −(2 − h)(3 − 5h + h^{2}),
C(h) = +6 − 19h + 17h^{2} − 3h^{3} = (2 − 3h)(3 − 5h + h^{2}), and
D(h) = −2 + 11h − 17h^{2} + 6h^{3} = −(1 − 2h)(1 − 3h)(2 − h).

Proof. Define h := n^{−1}, so h → 0 as n → ∞. From Corollary 9, there exists an asymptotic series of the form

−2σ_{n} ∼ ∑_{j≥0} d_{j}n^{−j−1/2}

as n → ∞. Moreover, d_{0} = 1. Define A(h) := (1 − h)(2 − 3h) in addition to B(h), C(h) and D(h). Using the recurrence (25) and the elementary identity 1/(n − m) = h/(1 − mh) for m ∈ {0, 1, 2, 3}, we have

∑_{j≥0} d_{j}(A(h)h^{j+1/2} + B(h) (\frac{h}{1 − h})^{j+1/2}

+ C(h) (\frac{h}{1 − 2h})^{j+1/2} + D(h) (\frac{h}{1 − 3h})^{j+1/2} ) ∼ 0.

Now, dividing both sides by h^{1/2}, we obtain

∑_{j≥0} d_{j}h^{j}(A(h) + B(h)(1 − h)^{−(j+1/2)}

+ C(h)(1 − 2h)^{−(j+1/2)} + D(h)(1 − 3h)^{−(j+1/2)}) ∼ 0.  (28)

An easy computation shows that

A(h) + B(h) + C(h) + D(h) = −4h^{2} + O(h^{3}),
B(h) + 2C(h) + 3D(h) = 8h + O(h^{2}), and
B(h) + 2^{2}C(h) + 3^{2}D(h) = O(h).

Thus, for all j ≥ 1, the terms involving d_{j} in (28) are 8jh^{j+2}+O(h^{j+3}). (The “8j” arises from −4 + 8(j + 1/2) = 8j.) This shows that the choice of d_{k} in (27) is necessary and sufficient to give an asymptotic series of the required form. Finally, we note that [h^{k+2−j}]A(h) = 0, since j ≤ k − 1 and deg(A(h)) = 2. Thus, a term involving A(h) has been omitted from (27).

Using Lemma 14, we computed the sequences (d_{n}) and (r_{n}) for n ≤ 1000, and verified the values previously computed (more slowly) via Corollary 9.

Since the power series occurring in (27) have a simple form, we can extract the coefficients of the required powers of h to obtain a recurrence for the d_{k}, as in Corollary 15. This gives a third way to compute the sequence (d_{n}).

Corollary 15. We have d_{0} = 1 and, for all k ≥ 1,

8k d_{k} = ∑_{j=0}^{k−1} α_{j,k} d_{j}.

Here

α_{j,k} = (−1 + 3 · 2^{m−1} − 2 · 3^{m})(τ)_{m−1}/(m − 1)!
+ (7 − 17 · 2^{m} + 17 · 3^{m})(τ)_{m}/m!
+ (−13 + 38 · 2^{m} − 33 · 3^{m})(τ)_{m+1}/(m + 1)!
+ 6(1 − 4 · 2^{m} + 3 · 3^{m})(τ)_{m+2}/(m + 2)!,  (29)

where m := k − j and τ := j + 1/2.

Proof (sketch). To prove Corollary 15, we apply the binomial theorem to the power series in (27), multiply by the polynomials B(h), C(h), and D(h), and extract the coefficient of h^{k+2−j}.

The following corollary is an easy deduction from Corollary 15, and gives an explicit recurrence for r_{k} = 2^{6k}d_{k}.

Corollary 16. We have r_{0} = 1 and, for all k ≥ 1,

k r_{k} = ∑_{j=0}^{k−1} β_{j,k} r_{j}, where β_{j,k} = 8^{2k−2j−1} α_{j,k} .

Although we have not proved Conjecture 10, the following result goes part of the way.

Theorem 17. For all k ≥ 0, we have k! r_{k} ∈ ℤ.

Proof. Let R_{k} := k!r_{k}. We show that R_{k} ∈ ℤ. From Corollary 16, R_{0} = 1 and, for k ≥ 1, R_{k} satisfies the recurrence

R_{k} = ∑_{j=0}^{k−1} β_{j,k} R_{j} \frac{(k − 1)!}{j!} .  (30)

The ratio of factorials in (30) is an integer, since j ≤ k − 1. Thus, in order to prove the result by induction on k, it is sufficient to show that β_{j,k} ∈ ℤ. Now, elementary number theory shows that 4^{ℓ}(j + 1/2)_{ℓ}/ℓ! ∈ ℤ for all j, ℓ ≥ 0. Thus, the expressions of the form (τ)_{m+δ}/(m + δ)! in (29) are in ℤ provided that m + δ ≥ 0. This is true as m ≥ k − j ≥ 1 and δ ≥ −1. To show that β_{j,k} ∈ ℤ, it is sufficient to have 8^{2m−1} ≥ 4^{m+2}, which holds for all m ≥ 2. In the case m = 1, it is easy to see that all the terms in (29) are in ℤ/4, so β_{m−1,k} = 8α_{m−1,k} ∈ ℤ. Thus, β_{j,k} ∈ ℤ for 0 ≤ j < k, and the result follows by induction on k.

Remark 18. The proof actually shows that β_{j,k} ∈ 2ℤ, which implies that R_{k} ∈ 2ℤ for all k > 0.

7 Acknowledgments

We thank Bruno Salvy for communicating his conjecture to one of us. An anonymous referee made helpful suggestions regarding the exposition. The first author was supported in part by ARC grant DP140101417.

References

[1] A. I. Aptekarev, On linear forms containing the Euler constant, preprint, 2009, http://

    arxiv.org/abs/0902.1768v2.

[2] P. Bala, Continued fraction for 1 − G, entry in OEIS \underline{A073003}, Oct. 11, 2012.

[3] R. P. Brent, A. J. Guttmann, and M. L. Glasser, A conjectured integer sequence arising from the exponential integral, preprint, 2018, http://arxiv.org/abs/1812.0316v2.

[4] R. P. Brent and F. Johannson, A bound for the error term in the Brent-McMillan algorithm, Math. Comp. 84 (2015), 2351–2359.

[5] A. R. Conway, A. J. Guttmann, and P. Zinn-Justin, 1324-avoiding permutations revisited, Adv. in Appl. Math. 96 (2018), 312–333.

[6] S. Covo, comment in OEIS \underline{A000262}, Jan. 20, 2010.

[7] A. Erdélyi (editor), Higher Transcendental Functions, Vol. 1, California Institute of Technology Bateman Manuscript Project, McGraw-Hill, 1953.

[8] L. Fejér, Sur une méthode de M. Darboux, C. R. Acad. Sci. Paris 147 (1908), 1040– 1042.

[9] P. Flajolet and R. Sedgewick, Analytic Combinatorics, Cambridge Univ. Press, 2009.

[10] G. Frobenius, Ueber die integration der linearen Differentialgleichungen durch Reihen, J. reine angew. Math. 76 (1873), 214–235.

[11] W. Gautschi, Computational aspects of three-term recurrence relations, SIAM Review 9 (1967), 24–82.

[12] A. J. Guttmann, Analysis of series expansions for non-algebraic singularities, J. Phys. A: Math. Theor. 48 (2015), 045209 (33pp).

[13] E. E. Kummer, De integralibus quibusdam definitis et seriebus infinitis, J. reine angew. Math. 17 (1837), 228-242.

[14] J. C. Lagarias, Euler’s constant: Euler’s work and modern developments, Bull. Amer. Math. Soc. 50 (2013), 527–628.

[15] National Institute of Standards and Technology, NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, version 1.0.19, 22 June 2018.

[16] A. Owczarek, T. Prellberg, and R. Brak, New scaling form for the collapsed polymer phase, Phys. Rev. Lett. 71 (1993), 951–953.

[17] O. Perron, Über das infinitäre Verhalten der Koeffizienten einer gewissen Potenzreihe, Archiv d. Math. u. Phys. (3), 22 (1914), 329–340.

[18] T. Rivoal, On the arithmetical nature of the values of the Gamma function, Euler’s constant, and Gompertz’s constant, Michigan Math. J. 61 (2012), 239–254.

[19] R. P. Stanley, Differentiably finite power series, Europ. J. Combinatorics 1 (1980), 175–188.

[20] N. M. Temme, Numerical evaluation of functions arising from transformations of formal series, J. Math. Anal. Appl. 51 (1975), 678–694.

[21] N. M. Temme, Remarks on Slater’s asymptotic expansions of Kummer functions for large values of the a-parameter, Adv. Dyn. Syst. Appl. 8 (2013), 365–377.

[22] F. G. Tricomi, Funzioni ipergeometriche confluenti, Edizioni Cremonese, Rome, 1954.

[23] D. N. Tulyakov, A system of recurrence relations for rational approximations of the Euler constant, Mat. Zametki 85 (2009), 782–787 (in Russian). English translation in Math. Notes 85 (2009), 746–750.

[24] M. Wallner, A bijection of plane increasing trees with relaxed binary trees of right height at most one, Theoret. Comput. Sci. 755 (2019), 1–12.

[25] E. M. Wright, The coefficients of a certain power series, J. Lond. Math. Soc. 7 (1932), 256–262.

[26] M. Wyman, The asymptotic behaviour of the Laurent coefficients, Canad. J. Math. 11 (1959), 534–555.

2010 Mathematics Subject Classification: Primary 34E05; Secondary 11Y55, 33C10, 33C15, 33F99. Keywords: asymptotics, confluent hypergeometric function, D-finite, Euler-Gompertz constant, exponential integral, Hadamard product, holonomic, Kummer function, modified Bessel function, stretched exponential.

(Concerned with sequences \underline{A000012}, \underline{A000262}, \underline{A067653}, \underline{A067764}, \underline{A073003}, \underline{A201203}, \underline{A293125}, \underline{A321937}, \underline{A321938}, \underline{A321939}, \underline{A321940}, \underline{A321941}, and \underline{A321942}.)

Received December 14 2018; revised versions received December 15 2018; May 27 2019. Published in Journal of Integer Sequences, August 15 2019.

Return to Journal of Integer Sequences home page.