Is the sequence of Apéry numbers a Stieltjes moment sequence?

Question

Consider the sequence of Apéry numbers
$$ A_n = \sum_{k=0}^n \binom{n}{k}\binom{n+k}{k}\sum_{j=0}^k \binom{k}{j}^3 = \sum_{k=0}^n \binom{n}{k}^2\binom{n+k}{k}^2 . $$ In an email, physicist Alan Sokal conjectures that it is a Stieltjes moment sequence. That is, that there exists a probability measure $\mu$ on $[0,+\infty)$ so that $$ A_n = \int_{[0,+\infty)} s^n\;d\mu(s) \tag{1}$$ for $n = 0,1,2,\dots$. [Of course you can equivalently say that $\mu$ is a nondecreasing function with $\mu(0)=0$ and $\lim_{x\to+\infty} \mu(s) = 1$ and that (1) is a Stieltjes integral, rather than a "measure" integral.]

Is that conjecture correct? Is $A_n$ a Stieltjes moment sequence?

[This question is a follow-up to A conjectured formula for Apéry numbers , where a formula for $A_n$ was established.]

Answer 1

Here are some computations that may be useful. Here below, $\big[ \,f\,\big]_a^b$ denotes $f(b)-f(a)$, and $n^{(k)}:=n(n-1)\dots(n-k+1)$.

1. (An integration by parts). Let $w$ be a solution of the third order linear ODE on the interval $[a,b]$: $$(x^4-34x^3+x^2)w''' + 3(2x^3-51x^2+x)w''+(7x^2-112x+1)w'+(x-5)w=0,$$ and put $$M(n):=\int_a^b x^nw(x)dx$$ for any integer $n\ge0$. Then, for any $n\ge 2$, $$n^3M(n)- (34n^3-51n^2+27n-5)M(n-1) + (n-1)^3M(n-2)= \Big[A(x)(x^n)'' + B(x)(x^n)' + C(x)x^n \Big]_a^b,$$ where $$ A:=pw \qquad B:= qw- (pw)' \qquad C:= (pw )'' - (qw)'+ rw ,$$ and $$p(x):=x^3-34x^2+x\qquad q(x):=3x^2-51x\qquad r(x):=x+10-x^{-1}.$$

Proof. We express the above polynomials of $n$ in the base $n^{(k)}$, then we absorb the latter terms as coefficients of derivatives of $x^n$, and finally we integrate by parts.

We have: $$n^3M(n)- (34n^3-51n^2+27n-5)M(n-1) + (n-1)^3M(n-2)=$$ $$\int_a^b\Big\{n^3x^nw -(34n^3-51n^2+27n-5)x^{n-1}w +(n-1)^3x^{n-2}w\Big\}dx=$$ $$\int_a^b\Big\{(n^ {(3)}+3n^ {(2)}+ n)x^nw -(34n^{(3)}+51n ^{(2)}+10n-5)x^{n-1}w +(n^{(3)}+n-1)x^{n-2}w\Big\}dx=$$ $$\int_a^b\Big\{(x^n)'''pw+(x^n)''qw+(x^n)'rw+x^n(5x^{-1} -x^{-2} )w\}dx=$$ $$ \int_a^bx^n \Big\{-(pw)'''+(qw)''-(rw)'+ (5x^{-1} -x^{-2} )w\}dx +\Big[ (x^n)''pw- (x^n)'(pw)'+ x^n(pw)'' +(x^n)'qw-x^n(qw)' +rw \Big]_a^b=$$ $$ -\int_a^bx^{n-1} \Big\{(x^4-34x^3+x^2)w''' + 3(2x^3-51x^2+x)w''+(7x^2-112x+1)w'+(x-5)w\Big\}dx+$$ $$+\Big[ (x^n)''pw+ (x^n)'\big(qw- (pw)'\big)+ x^n\big( (pw )'' - (qw)'+ rw \big) \Big]_a^b=$$ $$=\Big[A(x)(x^n)'' + B(x)(x^n)' + C(x)x^n \Big]_a^b.\qquad\square$$

2.(Consequence). Let $w$ a solution of the above linear equation on $(0,c)\setminus\{c_0\}$, with $\int_0^c w(x)dx=1$, and and assume it verifies the following linear boundary conditions, expressed in terms of the above coefficients $A,B,C$:

i) $A(x)=o(1),\quad B(x)=o(1),\quad C(x)=O(1)$, as $x\to0$;

ii) $A(x),\quad B(x),\quad C(x)$, are continuous at $x=c_0$ ,

iii) $A(x)=o(1) ,\quad B(x)=o(1) ,\quad C(x)=o(1)$, as $x\to c$.

Then the corresponding $M(n)$ are the Apéry sequence.

Indeed, computing the integral on $[0,c]$ for $M(n)$ as limit of integrals on $[\epsilon, c_0-\epsilon]\cup[c_0+\epsilon, c-\epsilon]$ as $\epsilon\to0$, and applying the above integration by parts formula, one gets that $M(n)$ satisfy the Apéry's recurrence, with $M(0)=1$ (note that $M(1)=5M(0)$ follows from the recurrence as well).

rmk. This also include the case where $w$ vanishes identically on $(0,c_0)$, and the condition is simply that $A, B, C$ should vanish both at $x=c_0$ and at $x=c$.

3. (Positive solutions of the third order ODE). Assume $u$ solves the second order linear ODE $$(x^3-34x^2+x)u''+(2x^2-51x+1)u'+\frac{1}{4}(x-10)u=0.$$ Then $w:=u^2$ solves $$(x^4-34x^3+x^2)w''' + 3(2x^3-51x^2+x)w''+(7x^2-112x+1)w'+(x-5)w=0.$$

Proof. Put $$P:=x^4-34x^3+x^2\qquad Q:=\frac{x^2}{2}-5x,$$ so the equation for $u$ (multiplied by $2x$) writes: $$2Pu''+P'u'+Qu=0,$$ and we have $$0=(2Pu''+P'u'+Qu)'u+ 3(2Pu''+P'u'+Qu)u'=$$ $$=(2Pu'''+2P'u''+P'u''+P''u'+Qu'+Q'u)u+3(2Pu''+P'u'+Qu)u'=$$ $$=P(2u'''u+6u''u')+3P'(u''u+u'^2)+(P''+4Q)u'u+Q'u^2=$$ $$=Pw''' +\frac{3}{2}P'w''+\Big(\frac{P''}{2}+2Q\Big)w'+Q'w,$$ which is the above third order equation for $w$.

Answer 2

I think the conjecture is true.

Below is an outline of a proof strategy---however, carefully verification of the details remains. If I find more time, I can try to fill those in (or maybe someone else provides a different proof before that).

For $\{A_n\}$ to be a Stieltjes-Moment sequence, two matrices $\Delta$ and $\Delta'$ must be positive definite.

The matrix $\Delta$ is defined as \begin{equation*} \Delta := \begin{pmatrix} A_0 & A_1 & \cdots & A_n\\ A_1 & A_2 & \cdots & A_{n+1}\\ & \vdots & & \\ A_n & A_{n+1} & \cdots & A_{2n} \end{pmatrix},\quad\text{i.e.}\quad \Delta_{ij} = A_{i+j-2}, 1\le i,j \le n+1, \end{equation*} while the matrix $\Delta' := [\Delta'_{ij}] = [A_{i+j-1}]$ for $1 \le i,j \le n+1$.

We prove below that $\Delta$ is symmetric positive semidefinite (a brief additional argument should establish strict positivity, which is what is needed to ensure infinite support).

First, we write $A_n$ using slightly different notation: \begin{equation*} A_n = \sum_{k=0}^n a_{n,k}^2,\qquad a_{n,k} := \binom{n}{k}\binom{n+k}{k}. \end{equation*} Next, define the order-0 Schmidt numbers \begin{equation*} S_n := \sum_{k=0}^n a_{n,k}, \end{equation*} and consider the matrix $M$ formed like $\Delta$ except that instead of $A_n$ we use $S_n$. We begin by proving that $S_n$ is positive definite, in particular by showing that \begin{equation*} S_{i+j-2} = \langle \phi(i), \phi(j) \rangle, \end{equation*} for some $\phi$. A similar (though more involved) technique can be followed for $A_n$ (though, if we actually could represent $a_{i+j-2,k}$ as an inner product, then the proof for $A_n$ would follow immediately using the Schur-product theorem).

The key trick is to use the ``symmetric'' form of the Vandermonde-Chu identity: \begin{equation*} \binom{r+s}{k} = \sum_{p,q \ge 0; p+q=k}\binom{r}{p}\binom{s}{q}. \end{equation*} Applying this identity, we have \begin{eqnarray*} \binom{i+j-2}{k} &=& \sum_{p,q\ge 0, p+q=k} \binom{i-1}{p}\binom{j-1}{q}\\ \binom{i+j-2+k}{k} &=& \sum_{p,q\ge 0, p+q=k} \binom{i-1+k/2}{p}\binom{j-1+k/2}{q}. \end{eqnarray*} Since $\binom{n}{j} = 0$ for $j > n$, we drop the summation indices (unless needed), and obtain \begin{eqnarray*} M_{ij} &=& \sum_{k}\biggl( \sum_{\substack{p,q \ge 0\\ p+q=k}} \binom{i-1}{p}\binom{j-1}{q} \biggr) \biggl( \sum_{p,q \ge 0, p+q=k} \binom{i-1+k/2}{p}\binom{j-1+k/2}{q}\biggr)\\ &=& \sum_k\sum_{\substack{p,q\ge 0, p+q=k\\ r,s \ge0, r+s=k}}\binom{i-1}{p}\binom{i-1+k/2}{r}\binom{j-1}{q}\binom{j-1+k/2}{s}\\ &=& \sum_{p,r,q, s \ge 0}\psi(i; p,r) \psi(j; q,s),\\ &=& \langle \psi(i), \psi(j) \rangle, \end{eqnarray*} for suitably defined $\psi(i; p, r)$. This proves that $M$ is a Gram matrix, hence positive definite.

In a similar way, we can prove that $\Delta_{ij} = \sum_k a_{i+j-2,k}^2 = \langle \phi(i), \phi(j)\rangle$ for a suitable mapping $\phi$, thus establishing positive definiteness of $\Delta$.

Continuing along this path, we can similarly prove $\Delta'$ is also positive definite, which will then finally establish the conjecture.

Answer 3

Some information on Pietro's ODE: $$ (4x^3-136x^2+4x)u''+(8x^2-204x+4)u'+(x-10)u=0 \tag{1}$$ I will use this notation: $$ c := (1+\sqrt{2}\;)^4 = 17+12\sqrt{2} \approx 33.97056 , \\ c_o := \frac{1}{c} = 34-c = 17-12\sqrt{2} \approx 0.02944 , \\ a := 1-c^2 = -576-408\sqrt{2} \approx -1159.9991 , \\ q := -\frac{11317}{4}-234\sqrt{2} \approx -660.176 , \\ \alpha := \frac{3}{2}, \beta := \frac{3}{2}, \gamma := \frac{3}{2}, \delta := 1, \epsilon := \frac{3}{2} . $$ Maple converts $(1)$ to a Heun differential equation, and evaluates it in terms of the Heun functions. See DLMF for information on that. I will follow their notation. In interval $(c_o,c)$, two linearly independent solutions of $(1)$ are $$ u_1(x) := (x-c_o)^{1/2}(c-x)^{1/2} Hl\big(a,q;\alpha,\beta,\gamma,\delta;1-cx\big) \\ = \sqrt {-{x}^{2}+34\,x-1}\;{\it Hl} \left( -408\,\sqrt {2}-576,-234\, \sqrt {2}-{\frac {1317}{4}};\frac{3}{2},\frac{3}{2},\frac{3}{2},1,-17\;x-12\,x\sqrt {2}+1 \right) \\ u_2(x) := (c-x)^{1/2} Hl\big(a,(a\delta+\epsilon)(1-\gamma)+q;\alpha+1-\gamma,\beta+1-\gamma,2-\gamma,\delta;1-cx\big) \\ = \sqrt {-x+17+12\,\sqrt {2}}\; {\it Hl} \left( -408\,\sqrt {2}-576,-42-30\,\sqrt {2},1,1,\frac{1}{2},1,-17 \,x-12\,x\sqrt {2}+1 \right) $$ The endpoints $x=c_o$ and $x=c$ correspond to $1-cx$ at the Heun singularities $0$ and $a$, respectively.

Here is a graph of $u_1(x)$

At the left end, $x=c_o$, it goes to zero, but has a vertical tangent (like a square-root)

At the right end, $x=c$, it goes to a finite limit, but has a vertical tangent

The wiggles are not real, but show Maple's difficulty in evaluating close to the singularity.

Here is a graph of $u_2(x)$

At the left end $x=c_o$ it approaches a definite value $2^{7/4}3^{1/2}$ with a non-vertical tangent

and at the right end $x=c$ it approaches a definite value with vertical tangent

For some linear combinations of $u_1$ and $u_2$ the square-root term cancels and it approaches a definite limit at $x=c$ with non-zero tangent. But for those, the approach at $x=c_o$ has vertical tangent.

Answer 4

OK, here is my current thinking on Pietro's approach.

We use the differential equation $$ (x^4-34x^3+x^2)U'''(x) + 3(2x^3-51x^2+x)U''(x)+(7x^2-112x+1)U'(x)+(x-5)U(x)=0 \tag{ODE3}$$ Multiply by $x^k$, $k \ge 0$, then integrate by parts as much as possible. The result is $$ \int \!{x}^{k-1} \left( {k}^{3}{x}^{2}-34\,{k}^{3}x +3\,{k}^{2}{x}^{2}+{k}^{3}-51\,{k}^{2}x+3\,k{x}^{2}-27\,kx+{x}^{2}-5\, x \right) U \left( x \right) {dx} = \left( {x}^{k+2}{k}^{2}-34\,{x}^{1+k}{k}^{2}+{x}^{k+2}k+{x}^{k}{k}^{2 }-17\,{x}^{1+k}k+{x}^{k+2}-10\,{x}^{1+k} \right) U \left( x \right) + \left( -{x}^{1+k}k-{x}^{3+k}k+34\,{x}^{k+2}k+{x}^{1+k}+2\,{x}^{3+k}- 51\,{x}^{k+2} \right) U' \left( x \right) + \left( {x}^{k+4}-34\,{x}^{3+k}+{x}^{k+2} \right) U'' \left( x \right) \tag{A}$$ In particular, for $k=0$, $$ \int (x-5)U(x)\;dx = \left( {x}^{2}-10\,x \right) U \left( x \right) + \left( 2\,{x}^{3}- 51\,{x}^{2}+x \right) U' \left( x \right) + \left( {x}^{ 4}-34\,{x}^{3}+{x}^{2} \right) U'' \left( x \right) \tag{A0}$$ Write $Q_k(x)$ for the right-hand-side of (A) and $\int R_k(x)U(x)\;dx$ for the left side. We will be doing differences like $$ \big[Q_k(x)\big]_a^b := Q_k(b)-Q_k(a) $$ since that will equal the integral on the left $\int_a^b R_k(x)U(x)\;dx$.
We want to arrange a solution $U$ so that $\int_0^c R_k(x)U(x)\;dx=0$ for $k \ge 0$. As Pietro noted, this will give us the recurrence we want for the moments $M(k):=\int_0^c x^kU(x)\;dx$. And in particular from (A0) we would have $\int_0^c (x-5)U(x)\;dx = 0$ so that $M(1)=5M(0)$.

Now consider (ODE3). First look at solutions as we approach $x=c$ from the left. Three linearly independent solutions are asymptotically: $1$, $(c-x)^{1/2}$, $(c-x)$. Other terms are higher powers (integer and half-integer). My calculations show:

for the term $1$, we get $Q_k(c^-)=1/2\, \left( 24+17\,\sqrt {2} \right) \left( 24\,k+24+7\,\sqrt {2} \right) \left( 17+12\,\sqrt {2} \right) ^{k}$

for the term $(c-x)^{1/2}$, we get $Q_k(c^-)=0$, amazingly

for the term $(c-x)$, we get $Q_k(c^-)=\left( -9792-6924\,\sqrt {2} \right) \left( 17+12\,\sqrt {2} \right) ^{k}$

terms $(c-x)^2$ and higher all produce $Q_k(c^-)=0$.

Now look at solutions as we approach $x=c_o$ from the right. Three linearly independent soltuions are asymptotically: $1$, $(x-c_o)^{1/2}$, $(x-c_o)$.

for the term $1$, we get $Q_k(c_o^+)=1/2\, \left( 17\,\sqrt {2}-24 \right) \left( -24\,k-24+7\,\sqrt {2} \right) \left( 17-12\,\sqrt {2} \right) ^{k}$

for the term $(x-c_o)^{1/2}$, we get $Q_k(c_o^+)=0$

for the term $(x-c_o)$, we get $Q_k(c_o^+)=-12\,\sqrt {2} \left( 17-12\,\sqrt {2} \right) ^{k+2}$

for higher terms, we get $Q_k(c_o^+)=0$

Now look at solutions as we approach $x=c_o$ from the left. Three linearly independent soltuions are asymptotically: $1$, $(c_o-x)^{1/2}$, $(x-c_o)$.

for the term $1$, we get $Q_k(c_o^-)=1/2\, \left( 17\,\sqrt {2}-24 \right) \left( -24\,k-24+7\,\sqrt {2} \right) \left( 17-12\,\sqrt {2} \right) ^{k}$

for the term $(c_o-x)^{1/2}$, we get $Q_k(c_o^-)=0$

for the term $(x-c_o)$, we get $Q_k(c_o^-)=-12\,\sqrt {2} \left( 17-12\,\sqrt {2} \right) ^{k+2}$

for higher terms, we get $Q_k(c_o^-)=0$

Now look at terms as we approach $x=0$ from the right. Three linearly independent terms are asymptotically: $1$, $\log x$, $(\log x)^2$.

for the term $1$, we get $Q_k(0^+)=0$

for the term $\log x$, we get $Q_k(0^+)=0$

for the term $(\log x)^2$, we get $Q_k(0^+)=0$ for $k \ge 1$ but $Q_0(0^+)=2$

for higher terms (such as $x$, $x\log x$, etc.), we get $Q_k(0^+)=0$

After all that, we come to the strategy for a solution. Start with the solution of (ODE3) approaching $c$ from the left that is asymptotically $(c-x)^{1/2}$. If we use only that, then $Q_k(c^-)=0$. (At the end we will multiply by a constant to get our final result.) Follow this solution to the left to the next singularity, that is approaching $c_o$ from the right. It will look like $r_1 + r_2(x-c_o)^{1/2}+r_3(x-c_o)+\dots$ for some constants $r_1, r_2, r_3$. Now consider solutions approaching $0$ from the right. Using only the two basis elements asymptotically $1$ and $\log x$, following them to the right until we approach $c_o$ from the left, we can choose a linear combination so that we get $r_1$ and $r_3$ from before: $r_1+r_2^*(c_o-x)^{1/2}+r_3(x-c)+\dots$. We want the same $r_1$ and $r_3$ as from the right, but $r_2^*$ is probably not the same as $r_2$. This way we get $Q_k(c_o^+)=Q_k(c_o^-)$ and $Q_k(0^+)=0$. So that all boundary terms cancel as desired.

Now we still have a constant factor to set. If $\int_0^c U(x)\;dx$ turns out to be nonzero, choose that factor so that $\int_0^c U(x)\;dx = 1$. Then (as noted at the beginning) we will have $A_n = \int_0^c x^n U(x)\;dx$ for all $n \ge 0$. Another thing we have to hope for is that $U(x)$ turns out to have constant sign: The way we constructed it, we have $U(x)>0$ near $c$, but we have to hope it turns out to be nonnegative everywhere.

I will see if I can get Maple to produce some pictures. I'm doing fine from $c$ going left to $c_o$, and from $c_o$ going both ways, but Maple seems reluctant to do (ODE3) from $0$ going up.

added
I am no expert on numerical ODEs. But here is what I get. Starting at the right, using the solution of (ODE3) with $U(x) = (c-x)^{1/2}+O((c-x)^{3/2})$ as $x \to c^-$, when I reach $c_o$ it looks like $$ U(x) = 1831.7 - 5769.99(x-c_o)^{1/2}-13916.8 (x-c_o)+O((x-c_o)^{3/2}) \qquad\text{as } x \to c_o^+ $$ Then, starting with the soution satisfying $$ U(x) = 1831.7 + 0.000103 (c_o-x)^{1/2}-13916.8 (x-c_o)+O((c_o-x)^{3/2}) \qquad\text{as } x \to c_o^- $$ we end up with $0$ on the $(\log x)^2$ term: $$ U(x) = O(\;|\log x|\;) \qquad\text{as } x \to 0^+ $$ The graph seems to be positive everywhere! So we merely need to divide by the integral.

Here is a plot. (Not very informative.)

Here it is near the point $c_o$: