Is the sequence of Apéry numbers a Stieltjes moment sequence? 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.]
 A: 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.
A: 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.
