Enumeration of lattice paths of a specific type

It seems the first solution to this problem appeared in Theorem 4 of

Raschel, Kilian, Counting walks in a quadrant: a unified approach via boundary value problems, J. Eur. Math. Soc. (JEMS) 14, No. 3, 749-777 (2012). ZBL1238.05014.

Notice that the walk corresponds to the fourth singular walk of Figure 2. The kernel $K(x,y;z)$ of the walk is given by $$K(x,y;z) = x y z(x +x/y +y +y/x -1/z).$$ Let $X_0(y) = X_0(y;z)$ (resp. $Y_0(x)=Y_0(x;z)$) be one of the solutions $x$ (resp. $y$) to $K(x,y;z) = 0$. Then Theorem 4 states that $$Q(x,y;z) = \sum_{m,n\geq 1} x^{m-1}y^{n-1} \sum_{\ell\geq 0} z^l a_{m,n}(\ell) $$ satisfies $$ Q(x,y;z) = \frac{1}{K(x,y;z)}\left(z x^2Q(x,0;z)+z y^2 Q(0,y;z)-xy\right)$$ with $$ Q(x,0;z) = \frac{1}{zx^2} \sum_{p\geq 0} Y_0 \circ (X_0\circ Y_0)^{\circ p}(x;z)\,\left[ (X_0\circ Y_0)^{\circ p}(x;z) - (X_0\circ Y_0)^{\circ (p+1)}(x;z) \right].$$ Here $f^{\circ p}$ means $f\circ \cdots \circ f$ with $p$ occurrences of $f$. Notice that by symmetry $Q(0,y;z) = Q(y,0;z)$.

A quick check with Mathematica reproduces the table for $\ell=7$:

k = x y z (x + x/y + y + y/x - 1/z);
x0[y_] = x /. Solve[k == 0, x][[1]] // FullSimplify[#, y > 0 && x > 0] &;
y0[x_] = y /. Solve[k == 0, y][[1]] // FullSimplify[#, y > 0 && x > 0] &;
Series[SeriesCoefficient[1/k (Sum[y0[xp] (xp - x0[y0[xp]]) /. 
   xp -> Nest[x0[y0[#]] &, x, p], {p, 0, 2}] 
   // -x y + # + (# /. x -> y) &), {z, 0, 7}], {x, 0, 7}, {y, 0, 7}]

yields the output

$$ \begin{aligned} &\ \ \ \left(y+21 y^2+80 y^3+125 y^4+85 y^5+21 y^6+y^7+O\left(y^8\right)\right)\\+&x\ \left(1+28 y+139 y^2+254 y^3+210 y^4+76 y^5+7 y^6+O\left(y^8\right)\right)\\+&x^2 \left(21+139 y+306 y^2+308 y^3+140 y^4+21 y^5+O\left(y^8\right)\right)\\+&x^3 \left(80+254 y+308 y^2+168 y^3+35 y^4+O\left(y^8\right)\right)\\+&x^4 \left(125+210 y+140 y^2+35 y^3+O\left(y^8\right)\right)\\+&x^5 \left(85+76 y+21 y^2+O\left(y^8\right)\right)\\+&x^6 \left(21+7 y+O\left(y^8\right)\right)\\+&x^7+O\left(x^8\right) \end{aligned} $$

This is an extension of the accepted answer. Using it, and ideas from some other papers, I managed to obtain an alternative expression for the generating function which I find interesting.

As explained in that answer, the problem reduces to finding $Q(x,0,z):=\sum_{i,\ell}a_{n,1}(\ell)x^nz^\ell$.

Let $z=\frac1{q+q^{-1}}$ and $x=-\frac{(1-q) (1-q^2) t}{(1-t) (1-q t)q}$. Then, $$ x^2zQ(x,0,z)=-x+\frac{1-4z^2}{z^2}\sum _{k=1}^\infty \frac{(-1)^kq^kt}{\left(1-q^kt\right)^2}. $$

The series $\sum _{k=1}^\infty \frac{(-1)^kq^kt}{\left(1-q^kt\right)^2}$ seems to be related to elliptic functions, Jacobi forms and such stuff, but I cannot figure out how exactly. Maybe I will ask a separate question about it.

Few words about the proof. The expression for the generating function in the other answer can be equivalently described as follows. Suppose $x$, $y$ satisfy $(x^2+y^2+x^2y+xy^2)z=xy$. Then, $$ (*)\qquad\qquad\qquad\qquad\qquad\qquad x^2zQ(x,0,z)=xy-y^2zQ(y,0,z).\qquad\qquad\qquad\qquad\qquad\qquad $$ Now the curve $(x^2+y^2+x^2y+xy^2)z=xy$ admits rational parametrizations $x(t)$, $y(t)$; moreover it is possible to choose such parametrization that satisfies $y(t)=x(qt)$, with $z=\frac1{q+q^{-1}}$ and $x(t)=-\frac{(1-q) (1-q^2) t}{(1-t) (1-q t)q}$. Then, denoting $F(t,q)=\frac{x(t)^2}{q+q^{-1}}Q\left(x(t),0,\frac1{q+q^{-1}}\right)$, we get from $(*)$ above the following functional equation for $F$: $$ F(t,q)=\frac{(1-q)^2 \left(1-q^2\right)^2 t^2}{(1-t) (1-q t)^2 \left(1-q^2 t\right)q}-F(qt,q). $$ This leads more or less straightforwardly to the claimed equality.

Update

Expression in terms of "original" variables can be also derived from this. Although less pretty, it is, I believe, still more explicit than the one in terms of iterated radicals.

We have $$ x^2zQ(x,0,z)=-x+\frac{1-4z^2}{z^2}\sum _{k=1}^\infty\frac{(-1)^k}{A(x,z)T_k\left(\frac{1}{2 z}\right)+B(x,z)U_{k-1}\left(\frac{1}{2 z}\right)-2}, $$ where $T$ and $U$ are the Chebyshev polynomials, while $A(x,z)$ and $B(x,z)$ are quadratic irrationalities determined by $$ \begin{aligned} A+\frac{2z}{1+2z}B&=2\left(1-\frac{1-2 z}{x z}\right),\\ A^2-\frac{4z^2}{1-4z^2}B^2&=4 \end{aligned} $$