By searching the Inverse Symbolic Calculator, we appear to be able to make the following conjecture about a real root to the equation:

$$\sum\limits_{k=0}^d \frac{x^k a_{k+1}}{k!}=0 \tag{1}$$

Let the lower triangular matrix $A$ be:

$$A=\binom{n-1}{k-1} a_{n-k+1} \tag{2}$$

where $n=1,2,3,4,5,...N$ and $k=1,2,3,4,5,...N$, with $N>>d$, and where the parentheses is the binomial function.

Calculate the matrix inverse $$B=A^{-1} \tag{3}$$, consider the first column of matrix $B$:

$$b_n=B(n,1) \tag{4}$$

and take the limiting ratio:

$$x=\lim_{n\to \infty } \, \frac{(n-1) b_{n-1}}{b_n} \tag{5}$$

Under what conditions for the coefficients:

$$a_1,...,a_{d} \tag{6}$$

is the limiting ratio $x$ in $(5)$ a real root solution to $(1)$

$$\sum\limits_{k=0}^d \frac{x^k a_{k+1}}{k!}=0$$ ?

Is the conjecture true at all?

I apologize for not letting the index of $a$ begin with $0$ instead of $1$.

Here is the Mathematica program for the conjecture:

```
Clear[A, B, a, b, x];
a = {1, 3, 5, 8, 5, 41, 39, 57, 53, 47, 13, 19, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0};
nn = Length[a];
d = Max[Flatten[Position[Sign[Abs[a]], 1]]]
A = Table[
Table[If[n >= k, Binomial[n - 1, k - 1]*a[[n - k + 1]], 0], {k, 1,
Length[a]}], {n, 1, Length[a]}];
b = Inverse[A][[All, 1]];
x = N[Table[(n - 1)*b[[n - 1]]/b[[n]], {n, nn - 8, nn - 1}], 30]
Sum[a[[k + 1]]/k!*x^k, {k, 0, d}]
Clear[x];
Sum[a[[k + 1]]/k!*x^k, {k, 0, d}]
```

The limiting ratio output from the program is: $$x=-0.474390307209018254579812222047$$

and this appears to be a solution to:

$$1+3 x+\frac{5 x^2}{2}+\frac{4 x^3}{3}+\frac{5 x^4}{24}+\frac{41 x^5}{120}+\frac{13 x^6}{240}+\frac{19 x^7}{1680}+\frac{53 x^8}{40320}+\frac{47 x^9}{362880}+\frac{13 x^{10}}{3628800}+\frac{19 x^{11}}{39916800}=0$$

In the program one needs to have $N$ much bigger than $d$ in order to see the conjecture in the output. Therefore there are a lot of trailing zeros in the coefficients $a$. You can add more zeros manually to the vector $a$ if you want to.

I don't know how to tag this question properly.

A related proof by joriki:

https://math.stackexchange.com/a/60385/8530

OEIS entry:

https://oeis.org/A167196

Related limiting ratio:

https://oeis.org/A132049

OEIS searches:

https://oeis.org/A322262

https://oeis.org/A006153

A much more conventional formulation of what I am doing in Mathematica:

```
Clear[x, b];
polynomial = (1 + 2*x + x^2/2! + x^3/3! + x^4/4! + x^5/5!);
digits = 100;
b = With[{nn = 200},
CoefficientList[Series[1/polynomial, {x, 0, nn}],
x] Range[0, nn]!] ;
nn = Length[b] - 10;
x = N[Table[(n - 1)*b[[n - 1]]/b[[n]], {n, nn - 8, nn - 1}], digits]
polynomial
```

Thanks to Harvey Dale in the OEIS for how to write the part:

```
b = With[{nn = 200},
CoefficientList[Series[1/polynomial, {x, 0, nn}],
x] Range[0, nn]!] ;
```

And thereby the polynomial need not have factorials:

```
Clear[x, b];
polynomial = (1 - 2 x + 3*x^2 - 5 x^3 + 7 x^4 - 11 x^5);
digits = 100;
b = With[{nn = 4000},
CoefficientList[Series[1/polynomial, {x, 0, nn}],
x] Range[0, nn]!] ;
nn = Length[b] - 10;
x = N[Table[(n - 1)*b[[n - 1]]/b[[n]], {n, nn - 8, nn - 1}], digits]
polynomial
```

The following program uses the method in the question by first Taylor expanding the Riemann zeta function at $zi$ and real part equal to $0$, and then adding trailing zeros to the vector $a$. Expanding at real part equal to $1$ gives a similar plot.

```
(*start*)
start = 10;
end = 30;
Monitor[list = Table[zz = 10;
d = 20;
a = Flatten[{CoefficientList[
Normal[Series[Zeta[x + z*N[I, d]], {x, 0, zz}]], x]*
Range[0, zz]!, Range[d]*0}];
nn = Length[a];
A = Table[
Table[If[n >= k, Binomial[n - 1, k - 1]*a[[n - k + 1]], 0], {k,
1, Length[a]}], {n, 1, Length[a]}];
Quiet[b = Inverse[A][[All, 1]]];
z*I + N[(nn - 1)*b[[nn - 1]]/b[[nn]], 40], {z, start, end, 1/10}],
z*10]
ListLinePlot[Re[list], PlotRange -> {-1, 3}, DataRange -> {start, end}]
ListLinePlot[Im[list], DataRange -> {start, end}]
(*end*)
```

The result is an approximation to the Riemann zeta zeros, where the plot of the real part stays around $\frac{1}{2}$ except at singularities at Gram points:

The heights of the steps in the staircase in the second plot are at imaginary parts of the Riemann zeta zeros.

The plots below are the same as above but from $z=10$ to $z=60$: