Extracting an asymptotic from a sequence defined by a recurrence relation

Suppose I have a sequence defined via its first term and a recurrence relation involving summation over all previous values with some coefficients. Here is the sequence I am interested in right now (although I would like to find a general method applicable to other sequences as well, if possible): $$a_0=1,\quad a_n=\sum_{k=0}^{n-1}\frac{a_k}{(n-k+1)!\,(2^n-1)}.\tag1$$ This sequence is positive rational, monotone decreasing, decaying faster than exponentially, with a few initial terms: $$1,\,\frac12,\,\frac5{36},\,\frac1{36},\,\frac{143}{32400},\,\frac{19}{32400},\,\frac{1153}{17146080},\,\frac{583}{857 30400},\,\dots\tag2$$

It is related to values of the Fabius function. I'm interested in its asymptotic expansion for $n\to\infty$. Using empirical numeric methods, and with some luck, I came to this conjecture: $$\log_2a_n\stackrel{\color{#aaaaaa}?}=-n \log_2n+\frac{n}{\log2}-\frac{\log _2^2n}{2} +\left(1-\frac{1}{2 \log ^22}-\frac{2}{\log2}\right)+\mathcal O\!\left(\frac{\log^2n}{n}\right)\!.\tag3$$ Could you suggest a way to prove it and find a few next terms of this asymptotic expansion?

How should I approach problems like this in general?

Update: The sequence in question is directly connected to the values of the Fabius function at negative integer powers of two: $$F\!\left(2^{-n}\right)=2^{-\binom n2}a_n.\tag4$$

There is also a direct non-recursive formula for it (conjectured): $$a_n\stackrel{\color{#aaaaaa}?}=\frac{(-1)^n}{(2;2)_n}\sum_{k=0}^n\frac{\binom n k_2}{2^{(n-1)k}(n+k)!}\sum_{\ell=0}^{2^k-1}(-1)^{\sigma_2(\ell)}\left(\ell-2^k+\tfrac12\right)^{n+k},\tag5$$ where $\left(q;q\right)_n=\prod_{k=1}^n\!\left(1-q^k\right)$ is the q‑Pochhammer symbol, ${\binom n k}_q=\frac{\left(q;q\right)_n}{\left(q;q\right)_k\,\left(q;q\right)_{n-k}}$ is the q‑binomial coefficient, and $\sigma_2(\ell)$ is the sum of binary digits of $\ell$; note that $(-1)^{\sigma_2\left(\ell\right)}$ is just the signed version of the Thue–Morse sequence.

Solution 1:

The following might help, but it's not a complete solution.

The coefficients $a_n$ have an alternative expression, where the $B_n$ are the Bernoulli numbers: $$ \tag{1} a_n = \sum_{k=0}^n 2^k a_k \ \frac{B_{n-k}}{(n-k)!} $$ To prove $(1)$, start with the recursion given in the problem statement, flip the direction of summation, and pull the $a_n$ term into the sum to get $$2^n a_n = \sum_{k=0}^n \frac{1}{k+1} \binom{n}{k} a_{n-k} \frac{(n-k)!}{n!} $$ Set $c_{n} = n!a_n$ so that $$2^n a_n = \sum_{k=0}^n \frac{1}{k+1} \binom{n}{k} c_{n-k} $$ This equation is invertible, e.g. Binomial Transform of a Particular Binomial Sum

Reverting back to $a_n$ gives the formula $(1)$. Define a generating function $f(x) =\sum_{n \ge 0} a_n x^n$. Then $(1)$ implies the functional equation

$$ \tag{2} \quad f(x) = f(2x) \frac{x}{e^x -1} $$

Iterate on $(2)$, and formally the solution is

$$ \tag{3} \quad f(x) = \prod_{k=0}^\infty \frac{ \exp{\big(x \ 2^{-k-1}\big)} - 1}{x \ 2^{-k-1} } $$

I've calculated many of the $a_n$ and checked that $(3)$ is satisfied numerically when the product is truncated and $x$ is sufficiently small.

$(3)$ looks like a formidable equation to do asymptotic analysis on. $(2)$ looks like a trick with a Mellin transform might be advantageous. The closest thing I know of that looks like the proposed asymptotic expansion has to do with Rice integrals, and instead of Bernoulli numbers, binomials appear in the sum. See example $4$ in $[1]$.

Edit 1/13/2022. Summary: the asymptotics of $a_n$ can be approximated by a saddle point analysis of the large $x$ asymptotics of $f(x)$ given above. Here is a short algorithm for calculation:

$$ \tag{4} x_0=-\frac{1}{\log{2}}W_{-1}\big( -\log{2} \cdot 2^{-n} \big) $$ $$ \tag{5} C=2^{-1/12}\exp{\Big(\big( \gamma_1 + \gamma^2/2-\pi^2/12 \big) }/\log{2} \Big) $$ $$ \tag{6} F(x)=C\sqrt{x}\exp{\big(x -\frac{1}{2 \log{2}}\log^2 x \big)}$$ $$b(x)=x \frac{d}{dx}(x \frac{d}{dx} \log{(F(x)} ) $$ $$\tag{7} a_n^{A} \sim (2 \pi b(x_0))^{-1/2} F(x_0) x_0^{-n} $$ Notation: $W_{-1}(x)$ is the Lambert W function on the '-1 branch' and is implemented in Mathematica as ProductLog[-1,x]. The $\gamma$ is Euler's constant and $\gamma_1$ is the first Stieltjes constant, $\gamma_1 \approx$ -0.0728158.

Here is a table comparing the true $a_n,$ the approximation (above) $a_n^{A},$ and the proposed approximation, $a_n^{P} = 2^{ (eq.3)}$ where eq. 3 corresponds to that in the proposers write up.

$$ \begin{array}{c|lcr} n & a_n & a_n^{A} & a_n^{P} \quad \\ \hline 20 & 6.08905 \,\cdot\, 10^{-22} & 6.07191\,\cdot\, 10^{-22} & 9.3967\,\cdot\, 10^{-22} \\ 40 & 1.00137 \,\cdot\, 10^{-52}& 1.00047\,\cdot\, 10^{-52} & 1.3980\,\cdot\, 10^{-52} \\ 60 & 1.30647\,\cdot\, 10^{-87} & 1.30565\,\cdot\, 10^{-87} & 1.7222\,\cdot\, 10^{-87} \\ 80 & 3.13549\,\cdot\, 10^{-125} & 3.13398\,\cdot\, 10^{-125} & 3.9804\,\cdot\, 10^{-125} \end{array} $$ The approximation $a_n^{A}$ is excellent, about 3 significant digits for $n=20,$ and about 4 digits for $n=80.$ A comment before sketching the proof. The proposer's solution seems to have missed the constant, and otherwise the expansion seems to be an asymptotic expansion of the Lambert W function, $$ W(x) \sim \log{x} - \log{\log{x}} + \cdots $$

The first part of the proof is to show that, for large $x$ $$ \tag{4} \quad f(x) = \prod_{k=1}^\infty \frac{ \exp{\big(x \ 2^{-k}\big)} - 1}{x \ 2^{-k} } \sim F(x)\exp{\big(P(x) \big)} $$ where $F(x)$ is given by eqs. (5) & (6) and $P(x)$ is given by

$$P(x)=\sum_{k=-\infty, k \ne 0}^\infty\ \frac{1}{2\pi i k}\Gamma\big(1-\frac{2 \pi i k}{\log2} \big) \zeta\big(1-\frac{2 \pi i k}{\log2} \big) \exp{\big(x \, \frac{2 \pi i k}{\log2} \big)}.$$

The formulas presented are a slight modification of that presented in Example 9.3 in reference [2]. We are going to ignore $P(x),$ in subsequent analysis, mainly because the amplitudes of the oscillations are tiny, with magnitudes of $7.4\,\cdot\,10^{-7}, 4.4\,\cdot\,10^{-13}\cdots$ for $k=1,2...$ However, if you want even more precise asymptotics (3 or 4 digits is not enough?) then they have to be included.

Reference [2] deals with a problem that has the functional equation $$ g(x) = g(2x)\frac{x}{1-e^{-x}} $$ Note how close that is to the functional equation in (3). On taking logarithms we have $$\log{\big(f(x)\big)} = \sum_{k=1}^\infty \log{\Big( \frac{ \exp{\big(x \ 2^{-k}\big)} - 1}{x \ 2^{-k} } \Big) } $$ $$ = \sum_{k=1}^\infty \log{\Big(\exp{\big(x \ 2^{-k}\big)}\Big) } + \log{\Big( \frac{ 1-\exp{\big(-x \ 2^{-k}\big)} }{x \ 2^{-k} } \Big) } $$ $$=x+ \sum_{k=1}^\infty \log{\Big( \frac{ 1-\exp{\big(-x \ 2^{-k}\big)} }{x \ 2^{-k} } \Big) }$$ The infinite sum is precisely that given in ref. [2] for $g(x)$, so its asymptotics can be copied verbatim, which I have done. Very briefly, the technique is to find the Mellin transform for the sum, then find the asymptotics by moving the contour of integration over the function's singularities. Later on, the reference extracts the asymptotics for the coefficients of $g(x)$ by using a process called dePoissionization, but that requires a growth restriction on the function. The addition of $e^x$ for my $f(x)$ destroys that possibility. However, traditional saddle point works. In that situation, one can use eq. (7) for a generic function $F(x).$ How does one determine the saddle point? Look at the Cauchy representation of $a_n,$ $$ a_n = \frac{1}{2 \pi i} \oint \frac{F(x)}{x^{n+1}} dx$$ $$ \rightarrow \frac{1}{2 \pi i} \oint \exp{\big( x - n\log{x} - \frac{1}{2 \log{2}}\log^2{x} \big)} \frac{dx}{\sqrt{x}} $$ The saddle point $x0$ is that in which the derivative of the argument of the exponential is zero,

$$ \frac{d}{dx} \big(x - n\log{x} - \frac{1}{2 \log{2}}\log^2{x} \big)\overset{!}{=} 0 $$ The solution to this equation for $x=x_0$ is given in eq. (4).

$[1]$ Mellin transform and asymptotics: Finite differences and Rice's integrals, P. Flajolet & R. Sedgewick, Ther. Comp. Sci 144 (1995) 101-124.

$[2]$ Average Case Analysis of Algorithms on Sequences W. Szpankowski, Wiley-Interscience, 2001.

Solution 2:

This is not a proper answer since it involves several heuristics. I will use the generating function $f(x)$ found by @skbmoore. For convenience, let $g(x)=f(2x)$. I would like to use Hayman's formula for the asymptotics of the Maclaurin coefficients of certain entire functions. The first non-trivial assumption is that $g(x)$ is admissible in the sense of Hayman, i.e., we can use his formula. Note that the Maclaurin coefficients of $g(x)$ are $2^n a_n$.

First, we estimate the growth of the logarithmic derivative. It can be shown that \begin{align*} a(r) :\!&= r\frac{{g'(r)}}{{g(r)}} = 2r + \sum\limits_{k = 0}^\infty {\left( {\frac{{r/2^k }}{{e^{r/2^k } - 1}} - 1} \right)} \\ & = 2r - \frac{{\log r}}{{\log 2}} - \frac{1}{2} - \sum\limits_{k = 1}^\infty {\frac{{2^k r}}{{e^{2^k r} - 1}}} +c(r), \end{align*} where the function $c(r)$ is bounded and staisfies $c(r)=c(2r)$. Numerics suggests that $|c(r)|<2\times 10^{-5}$. Hayman defines the sequence $r_n$ via $a(r_n ) = n$. Using the above, I found $$ r_n = \frac{n}{2} + \frac{{\log n}}{{2\log 2}} - \frac{1}{4} + o(1) $$ as $n\to +\infty$. We also have $b(r) := ra'(r) = 2r + \mathcal{O}(1)$ for large $r$. Now \begin{align*} \log g(r) & = \int_1^r {\frac{{a(t)}}{t}dt} + \log g(1) \\ & = 2r - \frac{{\log ^2 r}}{{2\log 2}} - \frac{{\log r}}{2} + C_0 - \sum\limits_{k = 1}^\infty \log (1 - e^{ - 2^k r} ) + d(r), \end{align*} with some real constant $C_0$ and with a bounded function $d(r)$ which staisfies $d(r)=d(2r)$. Therefore, by Hayman's formula, \begin{align*} \log (2^n a_n ) & = - n\log r_n + \log g(r_n ) - \frac{1}{2}\log (2\pi ) - \frac{1}{2}\log b(r_n ) + o(1) \\ & = - n\log \left( {\frac{n}{2}} \right) + n - \frac{{\log ^2 n}}{{2\log 2}} + C_1 +d(r_n)+ o(1) \end{align*} as $n\to +\infty$ with some real constant $C_1$. Equivalently, $$ \log _2 a_n = - n\log _2 n + \frac{n}{{\log 2}} - \frac{{\log _2^2 n}}{2} + C_2 +d(r_n)+ o(1) $$ as $n\to +\infty$ with some real constant $C_2$.