I am trying to understand the action of the second Macdonald operator $D_{q,t}^2$ on power sum polynomials $p_\lambda$, i.e., express $D_{q,t}^2 p_\lambda$ in terms of other $p_\mu$. After many steps of simplification, I reduce the problem to calculating the following expression: $$ B_n(u,v) = \frac{1}{4} \sum_{s=0}^{n-2} \sum_{r=0}^s \sum_{p=0}^s (-1)^{s+p} (t^2-1)^{n-2-s} (t-1)^s \binom{s}{r}$$ $$(s_{(r+p+u-s,v-r)} + s_{(r+p+v-s,u-r)} + s_{(p+u-r,r+v-s)} + s_{(p + v-r,u+r-s)}) e_{s-p}$$ Here $s_{a,b}$ denotes the Schur polynomial indexed by the integer vector $(a,b)$ and $e_a$ denotes the elementary symmetric polynomial of degree $a$. When $u=v=0$, the answer is pretty nice: $$ B_n(0,0) = \frac{(t^2-1)^{n-1} - (t-1)^{n-1}}{t^2 -t}$$.
I am wondering if someone knows how to deal with convolution product of Schur and elementary polynomials as stated. One thing I tried was to use the Jacobi-Trudy identity to write $s_{(a,b)} = h_a h_b - h_{a+1} h_{b-1}$, where $h_\lambda$ is the complete symmetric polynomial indexed by the partition $\lambda$.
Then use the relation between the generating functions of elementary and complete symmetric polynomials $E(-t) H(t) =1$. Unfortunately the expression above as $p$ ranges from $0$ to $s$ is only a partial convolution. Maybe I am missing something, but I am eager to learn from experts. The similar computation for $D_{q,t}^1$ can be found in the paper by Diaconis and Ram: stat.stanford.edu/~cgates/PERSI/papers/100726macdonaldpoly.pdf
