Strange polynomial analog of the Bell numbers 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}
 A: 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

Are the $\,P_n(σ)\,$ polynomials known?

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

We say a sequence R is Jacobi generated by the sequences
U and V if R are the coefficients of the series expansion
of the Jacobi continued fraction, recursively defined by
m = 1 - V(k)*x - U(k)*x^p/m, starting m = 1 and terminating
with 1/m, k iterating downwards from a given length to 1.
p is some integer (in the classic case p = 2).

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

Does the generating function $\sum_{n \geq 0} P_n(\sigma)
 z^n$ have a nice form?

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".
