This question addresses a hierarchy of linear recurrences which arise from an attempt to generalize the Nekrasov-Okounkov formula to the Young-Fibonacci setting. A related posting
extensions of the Nekrasov-Okounkov formula
asks how one might try to extend the Nekrasov-Okounkov formula by replacing the Plancherel measure on the Young lattice $\Bbb{Y}$ with another ergodic, central measure. In this discussion, I want to instead replace the Young lattice $\Bbb{Y}$ by the Young-Fibonacci lattice $\Bbb{YF}$ which comes equipped with its own Plancherel measure in virtue of being a $1$-differential poset. Allow me to briefly review some basics of the Young-Fibonacci lattice before I state the putative $\Bbb{YF}$-version of the Nekrasov-Okounkov partition function.
Young-Fibonacci Preliminaries: Recall that a fibonacci word $u$ is a word formed out of the alphabet $\{1,2\}$. As a set $\Bbb{YF}$ is the collection of a (finite) fibonacci words and $\Bbb{YF}_n$ will denote the set of fibonacci words $u \in \Bbb{YF}$ of length $|u|=n$ where $|u|:= a_1 + \cdots + a_k$ and where $u=a_k \cdots a_1$ is the parsing of $u$ into its digits $a_1, \dots, a_k \in \{1,2 \}$. The adjective Fibonacci reflects the fact that the cardinality of $\Bbb{YF}_n$ is the $n$-th Fibonacci number. I will skip defining the poset structure on $\Bbb{YF}$ and instead I point the readers to the Wikipedia page https://en.wikipedia.org/wiki/Young–Fibonacci_lattice. Suffice it to say that when endowed with an appropriate partial order $\unlhd$ the set $\Bbb{YF}$ becomes a ranked, modular (but not distributive), $1$-differential lattice. R. Stanley's $1$-differential property (see https://en.wikipedia.org/wiki/Differential_poset) is key here because it implies that the function $\mu^{(n)}_\mathrm{P}: \Bbb{YF}_n \longrightarrow \Bbb{R}_{>0}$ defined by
\begin{equation} \begin{array}{ll} \mu^{(n)}_\mathrm{P}(u) &\displaystyle := \ {1 \over {n!}} \, \dim^2(u) \quad \text{where} \\ \dim(u) &\displaystyle := \ \# \left\{ \begin{array}{l} \text{all saturated chains $(u_0 \lhd \cdots \lhd u_n)$ in $\Bbb{YF}$} \\ \text{starting with $u_0 = \emptyset$ and ending at $u_n =u$} \end{array} \right\} \end{array} \end{equation}
is a strictly positive probability distribution on $\Bbb{YF}_n$ for each $n \geq 0$. Furthermore, these distributions are coherent in the sense that the ratios
\begin{equation} \tilde{\mu}_\mathrm{P}(u \lhd v) \ := \ {\mu^{(n+1)}_\mathrm{P}(v) \over {\mu^{(n)}_\mathrm{P}(u)}} \end{equation}
restrict to a probability distribution on the set of covering relations $u \lhd v$ (i.e. edges in the Hasse diagram of $\Bbb{YF}$) for any fixed $u \in \Bbb{YF}_n$. We refer to $\mu^{(n)}_\mathrm{P}$ as the Plancherel measure for $\Bbb{YF}_n$. If $S:\Bbb{YF} \longrightarrow \Bbb{R}_{\geq 0}$ is some statistic let $\langle S \rangle_n$ denote its expectation value with respect to the Plancherel measure, i.e.
\begin{equation} \langle S \rangle_n \ := \ \sum_{|u|=n} \, {\dim^2(u) \over {n!}} \, S(u) \end{equation}
We may visualize a fibonacci word $u \in \Bbb{YF}$ using a profile of boxes akin to the way one depicts a partition by its Young diagram. The following example with $u = 12112211$ should illustrate the concept of a Young-Fibonacci diagram clearly. For emphasis each digit of the fibonacci word $u$ is written directly underneath the corresponding column of boxes:
\begin{equation} \begin{array}{cccccccc} & \Box & & & \Box & \Box & & \\ \Box & \Box & \Box & \Box & \Box & \Box & \Box & \Box \\ 1 & 2 & 1 & 1 & 2 & 2 & 1 & 1 \end{array} \end{equation}
A Fibonacci word $u$ will be synonymous with its Young-Fibonacci diagram and $\Box \in u$ will indicate membership of a box. The hook length $\mathrm{h}(\Box)$ of a box $\Box \in u$ is defined to be $1$ whenever it is in the top row; otherwise $\mathrm{h}(\Box)$ equals $1$ plus the total number of boxes directly above it and to its right. For example the hook lengths of the boxes of $u = 12112211$ are indicated in the tableaux below:
\begin{equation} \begin{array}{cccccccc} & \boxed{1 \ \ } & & & \boxed{1 \ \ } & \boxed{1 \ \ } & & \\ \boxed{11} & \boxed{10} & \boxed{8 \ \ } & \boxed{7 \ \ } & \boxed{6 \ \ } & \boxed{4 \ \ } & \boxed{2 \ \ } & \boxed{1 \ \ } \end{array} \end{equation}
These graphical conventions allows us to reformulate the value of the Plancherel measure as a product of hook-lenghts, i.e.
\begin{equation} \mu^{(n)}_\mathrm{P}(u) \ = \ \prod_{\Box \, \in \, u} \, {n! \over {\mathrm{h}^2(\Box)} } \end{equation}
This is a non-trivial observation made by R. Stanley in the course of his work examining differential posets.
The $\Bbb{YF}$-version of the Nekrasov-Okounkov partition function: For a fibonacci words $u \in \Bbb{YF}$ define a $t$-statistic $H_t(u) := \prod_{\Box \, \in \, u} \, \big(\mathrm{h}^2(\Box) - t \big)$ and the $\Bbb{YF}$-Nekrasov-Okounkov partition function as
\begin{equation} \begin{array}{ll} F(z;t) &\displaystyle = \ \sum_{n \geq 0} {z^n \over {n!}} \, \langle H_t \rangle_n \\ &\displaystyle = \ \sum_{n \geq 0} {z^n \over {n!}} \, \sum_{|u|=n} \, {\dim^2(u) \over {n!}} \, H_t(u) \end{array} \end{equation}
It will be convenient, when dealing with expansions into elementary symmetric polynomials, to make the change of variable $z \mapsto -z$ and consider $F^\vee(z;t) := F(-z;t)$ instead; the effect of this sign-change is to replace the statistic $H_t(u)$ by $H^\vee_t(u) := \prod_{\Box \, \in \, u} \, \big(t -\mathrm{h}^2(\Box) \big)$ in the definition of the partition function. After expanding into elementary symmetric polynomial $E_k$ we get
\begin{equation} H^\vee_t(u) \ = \ \sum_{k=1}^n \, (-t)^{n-k} \, E_k \big( \mathrm{h}^2(\Box) \big)_{\Box \, \in \, u} \end{equation}
and
\begin{equation} F^\vee(z;t) \ = \ \sum_{k \geq 0} \, (-t)^{n-k} \, \overbrace{\sum_{n \geq 0} \, {z^n \over {n!}} \, \langle E_k \rangle_n}^{F^\vee_k(z)} \end{equation}
which effectively reduces the problem of calculating $F^\vee(z;t)$ to the problem of evaluating the expectation values $\langle E_k \rangle_n$.
Evaluating expectation values: Fibonacci words $u \in \Bbb{YF}_n$ with $n \geq 2$ can be separated into two disjoint groups: Those of the form $u=1v$ for $v \in \Bbb{YF}_{n-1}$ and those of the form $u=2v$ for $v \in \Bbb{YF}_{n-2}$. Depending on whether the prefix of $u$ is $1$ or $2$ we can write down a recursive formula for the value of $E_k(u) := E_k \big( \mathrm{h}^2(\Box) \big)_{\Box \, \in \, u}$ by analyzing the hook length(s) of the box(es) in the left-most column, specifically:
\begin{equation} \begin{array}{lll} E_k(1v) &= E_k(v) + n^2E_{k-1}(v) &\text{if} \ |v| = n-1 \\ E_k(2v) &= E_k(v) + (n^2+1)E_{k-1}(v) + n^2E_{k-2}(v) &\text{if} \ |v| = n-2 \end{array} \end{equation}
Using the observation that $\dim(1v) = \dim(v)$ and $\dim(2v) = (|v| + 1)^2 \dim(v)$ we may conclude
\begin{equation} \langle E_k \rangle_n = \left\{ \begin{array}{l} \displaystyle {1 \over n} \langle E_k \rangle_{n-1} \ + \ {n-1 \over n} \langle E_k \rangle_{n-2} \\ \\ \displaystyle + \ n \langle E_{k-1} \rangle_{n-1} \ + \ {(n-1)(n^2+1) \over n} \langle E_{k-1} \rangle_{n-2} \ + \ n(n-1) \langle E_{k-2} \rangle_{n-2} \end{array} \right. \end{equation}
If we set $\sigma_k(n) := {1 \over {n!}} \, \langle E_k \rangle_n$ then the above recursion can be rewritten as:
\begin{equation} n^2\sigma_k(n) = \underbrace{\sigma_k(n-1) + \sigma_k(n-2)}_{\text{homogeneous part}} + \underbrace{n^2\sigma_{k-1}(n-1) +(n^2 +1)\sigma_{k-1}(n-2) + n^2\sigma_{k-2}(n-2)}_{\text{inductive heap of inhomogeneous junk}} \end{equation}
which can be converted, using the usual yoga of generating functions, into the following second order inhomogeneous ODE for $F^\vee_k(z) := \sum_{n \geq 0} \sigma_k(n) z^n$
\begin{equation} z^2 \, {d^2 \over {dz^2}} \, F^\vee_k \ + \ z \, {d \over {dz}} \, F^\vee_k \ + \ \big(z^2 + z \big) \, F^\vee_k \ = \ G_{\leq k}(z) \ + \ \big( \sigma_k(1) - \sigma_k(0) \big)z \end{equation}
where $G_{\leq k}(z)$ is the generating function associated to the heap of inhomogeneous junk which, by induction, will have been previously evaluated. The homogeneous ODE has two nice independent solutions $Y_1(z) = e^z$ and $Y_2(z)= e^z \int z^{-1} e^{-2z} dz$ whose Wronskian is just $W={z^{-1}}$. One starts the inductive engine beginning with $F^\vee_0(z) = e^z$. For $k=1$ its not hard to see that $\sigma_1(1)=1$ and $\sigma_1(0)=0$ while
\begin{equation} \begin{array}{ll} \displaystyle G_{\leq 1}(z) &\displaystyle = \ {z^2 \over {1-z}} \ + \ 2 \sum_{n \geq 2} \, n^2 z^n \\ &\displaystyle = {z^2 \over {1-z}} \ + \ 2z \, \Bigg( {1+z \over {(1-z)^3}} \, - \, 1 \Bigg) \end{array} \end{equation}
so the ODE for $F^\vee_1(z)$ becomes
\begin{equation} z^2 \, {d^2 \over {dz^2}} \, F^\vee_1 \ + \ z \, {d \over {dz}} \, F^\vee_1 \ + \ \big(z^2 + z \big) \, F^\vee_1 \ = \ z + G_{\leq 1}(z) \end{equation}
By variation of parameters, a particular inhomogeneous solution is
\begin{equation} \begin{array}{rl} \displaystyle Y_\mathrm{particular}(z) &\displaystyle = \ V_1(z) \cdot e^z \ + \ V_2(z) \cdot e^z \overbrace{\int {dz \over z} \, e^{-2z} }^{\gamma(z)} \\ \displaystyle V_1(z) &\displaystyle = \ -\int z \, e^z \, \gamma(z) \, \Big(z + G_{\leq 1}(z) \Big) \, dz \\ \displaystyle V_2(z) &\displaystyle = \ \ \ \ \ \int z \, e^z \, \Big(z + G_{\leq 1}(z) \Big) \, dz \end{array} \end{equation}
Question. Has the linear recurrence satisfied by $\sigma_k(n)$ or else the hierarchy of 2nd order inhomogeneous ODEs been studied?
thanks, ines.