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

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}

Here is my heuristics understanding on why the summation is so special. Consider the function

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

Let $a \in \Bbb{C} \setminus \{0,1,2,\cdots\}$ and $\epsilon > 0$ be sufficiently small. Putting the matter of rigor aside, computing the following contour integral

$$ \frac{1}{2\pi i} \oint_{|z| = \epsilon} \frac{f(z)}{z-a} \, dz $$

in two ways yields

$$ \sum_{n=0}^{\infty} \frac{1}{n!}\frac{1}{n-a} = -e f(a). \tag{1}$$

Moreover, by the direct computation we check that

$$ f(z-1) = zf(z) - 1. \tag{2}$$

Now what is so special about our sum is that, if we write $\omega = \frac{1}{2} + i\frac{\sqrt{3}}{2}$, then applying $\text{(1)}$ to the partial fraction decomposition

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

we find that

$$\sum_{n=0}^{\infty} \frac{1}{n!} \frac{1}{n^4+n^2+1} = -\frac{e}{2i\sqrt{3}} \left[ \bar{\omega}^2 (\omega f(\omega) - f(\omega-1)) - \omega^2 (\bar{\omega} f(\bar{\omega}) - f(\bar{\omega} - 1)) \right]. $$

Now applying $\text{(2)}$ gets rid of the occurrence of the nebulous function $f$ and simply we get

$$\sum_{n=0}^{\infty} \frac{1}{n!} \frac{1}{n^4+n^2+1} = \frac{e}{2}. $$

As we see, all this argument relies on the special structure of the partial fraction decomposition $\text{(3)}$.

EDIT. By translating $\text{(2)}$ in terms of the left-hand side of $\text{(1)}$, we can give a more elementary proof.

Let $a \in \Bbb{C}\setminus\{0,1,2,\cdots\}$. Then

\begin{align*} \sum_{n=0}^{\infty} \frac{1}{n!} \left( \frac{1}{n-(a-1)} - \frac{a}{n-a} \right) &= \sum_{n=0}^{\infty} \frac{1}{n!} \frac{1}{n-(a-1)} - \sum_{n=0}^{\infty} \frac{1}{n!} \frac{a}{n-a} \\ &= \sum_{n=1}^{\infty} \frac{1}{(n-1)!} \frac{1}{n-a} - \bigg( -1 + \sum_{n=1}^{\infty} \frac{1}{n!} \frac{a}{n-a} \bigg) \\ &= 1 + \sum_{n=1}^{\infty} \frac{1}{n!} \bigg( \frac{n}{n-a} - \frac{a}{n-a} \bigg) \\ &= e. \end{align*}