Let $\vec{x} = (x_0, x_1, x_2, \dots)$ and $\vec{y}=(y_1,y_2,y_3, \dots)$ be two systems of parameters/variables. The Motzkin polynomials $P_n(\vec{x},\vec{y})$ for $n \geq 0$ are defined by the following quadratic recursion
\begin{equation} P_n(\vec{x},\vec{y}) \ = \ x_0 P_{n-1}(\vec{x},\vec{y}) \, + \, \sum_{i \, = \, 0}^{n-2} \, y_1P_i(\vec{x}_+, \vec{y}_+) P_{n-i-1}(\vec{x},\vec{y}) \end{equation}
where $P_0(\vec{x},\vec{y}) := 1$ and $P_1(\vec{x},\vec{y}) := x_0$ are the initial polynomials and where $\vec{x}_+:= (x_1, x_2, x_3, \dots)$ and $\vec{y}_+:=(y_2,y_3,y_4, \dots)$ are the respective one-step shifts of the sequences $\vec{x}$ and $\vec{y}$. From an enumerative standpoint, the Motzkin polynomial $P_n(\vec{x},\vec{y})$ is the multi-variate generating function of all Motzkin paths with $n$ steps: Each horizontal step (at level $k$) is weighted $x_k$, each ascending step (at level $k$) is weighted $y_k$, and the overall weight of the path is the product of weights of all horizontal steps and ascents which are taken. I should add that the generating function $\sum_{n \geq 0} P_n(\vec{x},\vec{y})z^n$ can be be formally expressed as the following Jacobi continued fraction involving the parameters $\vec{x}$ and $\vec{y}$
\begin{equation} J_{\vec{x}, \vec{y}} \, (z) \ := \ {1 \over {1 - x_0z - {\displaystyle y_1z^2 \over {\displaystyle 1 - x_1 z - {\displaystyle y_2z^2 \over {\displaystyle 1 - x_2 z - {\displaystyle y_3z^2 \over {\ddots } } } } } } } } \end{equation}
I'm concerned with the following specialization. For integers $k \geq 1$ let $x_{k-1} = \sigma + k - 1$ and $y_k = \sigma + k -1$ where $\sigma$ is an indeterminate. Since $\sigma$ is the operative variable I shall, for brevity's sake, write $P_n(\sigma)$ instead of $P_n(\vec{x},\vec{y})$. Brute force calculations reveal that
\begin{equation} \begin{array}{l} P_0(\sigma) \, =1 \\ P_1(\sigma) \, = \sigma \\ P_2(\sigma) \, = \sigma^2 + \sigma\\ P_3(\sigma) \, = \sigma^3 + 3\sigma^2 + \sigma \\ P_4(\sigma) \, = \sigma^4 + 6\sigma^3 + 6\sigma^2 + 2\sigma \\ P_5(\sigma) \, = \sigma^5 + 10\sigma^4 + 20\sigma^3 + 16\sigma^2 + 5\sigma \\ P_6(\sigma) \, = \sigma^6 + 15\sigma^5 + 50\sigma^4 + 71\sigma^3 + 51\sigma^2 + 15\sigma \\ P_7(\sigma) \, = \sigma^7 + 21\sigma^6 + 105\sigma^5 + 231\sigma^4 + 281\sigma^3 + 186\sigma^2 + 52\sigma \end{array} \end{equation}
Evidence seems to suggest that $P_n(1) = B_n$ where $B_n$ is the $n$-th Bell number while $P_n(2)$ counts the number of irreducible set partitions of size $n$ (see sequence A074664 at the OEIS). Furthermore
\begin{equation} \begin{array} \big[\sigma] \, P_n(\sigma) &\displaystyle = \, B_{n-2} \ \text{for $n \geq 2$} \\ \big[\sigma^{n-2}\big] \, P_n(\sigma) &\displaystyle = \, {1 \over {12}} \, n^2(n^2-1) \ \text{for $n \geq 3$} \\ \big[\sigma^{n-1}\big] \, P_n(\sigma) &\displaystyle = \, \binom{n}{2} \ \text{for $n \geq 2$} \\ \end{array} \end{equation}
where $[\sigma^k] \, P_n(\sigma)$ denotes the coefficient of $\sigma^k$ occurring in $P_n(\sigma)$. For $0 \leq n \leq 3$ we have $P_n(\sigma) = T_n(\sigma)$ where $T_n(\sigma)$ is the Touchard polynomial but this coincidence ceases for $n \geq 4$.
Question: Are the $P_n(\sigma)$ polynomials known? Does the generating function $\sum_{n \geq 0} P_n(\sigma) z^n$ have a nice form?
thanks, jeanne.
Up-date: Rephrasing the calculation in Somos' post and also Following Qiaochu Yuan's suggestion to use the continued fraction expansion of the generating function $J(\sigma, z):= \sum_{n \geq 0} P_n(\sigma) z^n$ we get the functional equation
\begin{equation} J(\sigma, z) \ = \ {1 \over {1 - \sigma z - \sigma z^2 J(\sigma +1 ,\ z) }} \end{equation}
which we can re-write in terms of linear fractional transformations as
\begin{equation} \begin{array}{ll} \begin{pmatrix} 1 - \sigma z & -1 \\ \sigma z^2 & 0 \end{pmatrix} \cdot J(\sigma, z) &\displaystyle = \ {(1 - \sigma z)J(\sigma,z) \, - \, 1 \over{\sigma z^2 J(\sigma, z)}} \\ \\ &= \ J(\sigma +1 ,z) \end{array} \end{equation}
The term "nice form" might entail having a differential-recursive formula for the polynomials $P_n(\sigma)$ analogous to the Rodrigues formula for the Touchard polynomials, namely:
\begin{equation} T_{n+1}(x) \ = \ x \Big(1 \, + \, {d \over {dx}} \Big) T_n(x) \end{equation}
