A nice infinite series: ${\sum\limits_{n=1}^\infty}\frac{1}{n!(n^4+n^2+1)}=\frac{e}{2}-1$ - looking for a more general method

Question

while clearing out some stuff I found an old proof I wrote for an recreational style problem a while back, here's the sum:

${\sum\limits_{n=1}^\infty}\frac{1}{n!(n^4+n^2+1)}=\frac{e}{2}-1$

I'll show how I got to the limit in a moment, but because my argument was a bit long winded / ad-hoc (first thing that came to mind at the time), just wondering if this kind of inverse polynomial/factorial sum can be attacked using other techniques, is it part of a well known / classical family of these kind of sums / any useful references? Thanks in advance!

Proof: \begin{align} \displaystyle{\sum\limits_{n=1}^\infty}\frac{1}{n!(n^4+n^2+1)} &= \displaystyle{\sum\limits_{n=1}^\infty}\frac{1}{n!{(n^2+1)^2-n^2}}= \displaystyle{\sum\limits_{n=1}^\infty}\frac{1}{n!(n^2+1-n)(n^2+1+n)}\nonumber \\ &= \displaystyle{\sum\limits_{n=1}^\infty \frac{n}{n \cdot n!(n^2+1-n)(n^2+1+n)} }\nonumber\\ &= \displaystyle{\frac{1}{2}\sum\limits_{n=1}^\infty\frac{1}{n \cdot n!} \left(\frac{1}{n^2+1-n}-\frac{1}{n^2+1+n}\right)}\end{align} Now this almost telescopes, the factorial term being the issue, but if we slightly change the groupings of the summands (essentially starting the summation having moved one index along) we get something that we can further simplify by taking the telescoping part out of the parentheses and putting the factorial part in:

$ = \displaystyle{\frac{1}{2}\left[\frac{1}{1\cdot1!}\cdot\frac{1}{1}+\sum\limits_{n=2}^{\infty}\frac{1}{n^2+1-n}\left\{\frac{1}{n \cdot n!}-\frac{1}{(n-1) (n-1)!}\right\}\right]}$

We have taken out the 1/quadratic that was common to successive telescoping summands in the original summation, and taken the difference of the factorials instead since as we will see this simplifies nicely:

$ \begin{align}&= \displaystyle{\frac{1}{2}\left[1+\sum\limits_{n=2}^{\infty}\frac{1}{n^2+1-n}\left\{\frac{n-1-n^2}{n\cdot(n-1)\cdot n!}\right\}\right]}\nonumber \\ &= \displaystyle{\frac{1}{2}\left[1+\sum\limits_{n=2}^{\infty}\frac{1}{n\cdot(n-1)\cdot n!}\right]}\end{align}$

Using the same partial fraction trick we used in (1) to pseudo-telescope and regroup terms to get to (2) helps again:

\begin{align} &= \displaystyle{\frac{1}{2}\left[1+\sum\limits_{n=2}^{\infty}\frac{1}{n!}\left(\frac{1}{n-1}-\frac{1}{n}\right)\right]} =\displaystyle{\frac{1}{2}\left[1-\left\{\frac{1}{2!}+\sum\limits_{n=3}^{\infty}\frac{1}{(n-1)}\left(\frac{1}{n!}-\frac{1}{(n-1)!}\right)\right\}\right]}\nonumber \\ &=\displaystyle{\frac{1}{2}\left[1-\left\{\frac{1}{2}-\sum\limits_{n=3}^{\infty}\frac{1}{n!}\right\}\right]} =\displaystyle{\frac{1}{2}\left[\frac{1}{2}+\left\{\sum\limits_{n=0}^{\infty}\frac{1}{n!}\right\}-2 \frac{1}{2}\right]}\\ &=\displaystyle{\frac{1}{2}\left\{\sum\limits_{n=0}^{\infty}\frac{1}{n!}\right\}-1=\frac{e}{2}-1}\,\,\,\square \end{align}

Answer 1

In this answer, we try to explain what is so special about OP's sum. Consider the function

$$ g(z) = \sum_{n=0}^{\infty} \frac{1}{n!}\frac{1}{n-z}. $$

Then

\begin{align*} g(z-1) - zg(z) &= \sum_{n=0}^{\infty} \frac{1}{n!} \frac{1}{n-(z-1)} - \sum_{n=0}^{\infty} \frac{1}{n!} \frac{z}{n-z} \\ &= \sum_{n=1}^{\infty} \frac{1}{(n-1)!} \frac{1}{n-z} - \bigg( -1 + \sum_{n=1}^{\infty} \frac{1}{n!} \frac{z}{n-z} \bigg) \\ &= 1 + \sum_{n=1}^{\infty} \frac{1}{n!} \bigg( \frac{n}{n-z} - \frac{z}{n-z} \bigg) \\ &= e. \tag{1} \end{align*}

Now, let $\omega = e^{i\pi/3}$. The key observation is that $\frac{1}{z^4 + z^2 + 1}$ admits partial fraction decomposition of the form

$$ \frac{1}{z^4 + z^2 + 1} = \frac{1}{2i\sqrt{3}} \left( \frac{\omega^2}{z - (\bar{\omega}-1)} - \frac{\omega}{z-\bar{\omega}} + \frac{\bar{\omega}}{z-\omega} - \frac{\bar{\omega}^2}{z - (\omega-1)} \right). \tag{2} $$

Hence it follows that

\begin{align*} \sum_{n=0}^{\infty} \frac{1}{n!} \frac{1}{n^4+n^2+1} &= \frac{1}{2i\sqrt{3}} \left[ \omega^2 (g(\bar{\omega} - 1) - \bar{\omega} g(\bar{\omega})) - \bar{\omega}^2 (g(\omega-1) - \omega g(\omega)) \right] \\ &= \frac{1}{2i\sqrt{3}} ( \omega^2 e - \bar{\omega}^2 e ) \\ &= \frac{e}{2}. \end{align*}

As we see, the argument critically hinges on the very special form of the decomposition $(2)$ that is perfectly tailored for the functional identity $(1)$.

---

Just for fun. In my original argument, I introduced an auxiliary function $f(z)$ defined by

$$ f(z) = \sum_{n=0}^{\infty} \frac{1}{z(z-1)\cdots(z-n)}. $$

Although the new answer does not make use of this $f(z)$, we may still investigate $f(z)$ on its own. One interesting property is that $f(z)$ is related to $g(z)$ by the relation

$$ g(z) + ef(z) = 0. $$

I originally derived this using contour integration, but I believe this can also be established by showing that LHS has only removable singularities and is bounded.