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}
Define the monic polynomials $\,P_n(x)\,$ of degree $n$ in $x$ such that its ordinary generating function $\, f_x(z) := \sum_{n=0}^\infty P_n(x)z^n \,$ satisfies the equation
$$ f_x(z) = \frac1{1 - xz - xz^2\,f_{x+1}(z)} = 1 + (x)z + (x+x^2)z^2 + (x+3x^2+x^3)z^3 + \cdots.$$
This equation is equivalent to the recursion
$$ P_0(x) = 1,\quad P_{n+1}(x) = x\,P_n(x) + x\sum_{k=0}^{n-1} P_k(x+1)P_{n-1-k}(x). $$
The answer to your first question
is probably "no" but should be "yes". A lookup in the OEIS did not find any match for the triangular sequence of coefficients which is strong evidence that it is not yet a known sequence but not conclusive. I created A357438 for this sequence and so that may contribute to its reknown. However, Peter Luschny created in 2018 A321960 which is the array $\,\{P_n(k)\}_{n,k=0}^\infty\,$ read by antidiagonals with a definition equivalent to $\,f_x(z)$ and it links to A321964 which states
In general, given a sequence of polynomials with integer coefficients, the triangular sequence of coefficients may be listed in the OEIS and/or the array of values of the polynomials for nonnegative integers read by antidiagonals. In this particular case, now both are in the OIES
My answer to your second question
is that this is hard to tell since "nice form" is not a well defined term. I am still looking for something like that.
You may be interested in the closely related blog article by Qiaochu Yuan "Moments, Hankel determinants, orthogonal polynomials, Motzkin paths, and continued fractions".