Since the Laplace operator $$ \Delta=\frac{\partial^2}{\partial x_1^2}+\cdots+\frac{\partial^2}{\partial x_n^2} $$ preserves symmetric polynomials in $n$ variables, its action on the Schur polynomial $s_{\lambda} $ can again be expressed in terms of Schur polynomials $$ \Delta(s_\lambda) =\sum_{\mu} c_{\lambda \mu} s_\mu, $$ for some $\mu $ and $c_{\lambda \mu}.$
For the case $n=2, a>b+1$, I have empirically established the following explicit formula $$ \frac{1}{2} \Delta(s_{a,b})=\frac{b(b-1)}{2}s_{a,b-2}+b s_{a-1,b-1}+\left(\frac{a(a-1)}{2}+b+1 \right)s_{a-2,b} -\sum_{i=b+1}^{\left[\frac{a+b-2}{2}\right]} (a-2i+b-1) s_{a-i+b-2,i}. $$
I tried, considering that the action of $\Delta $ on the power sums $p_k$ is quite simple, to use the well-known decomposition of Schur polynomials through power sums (Murnaghan-Nakayama Rule) but so far without much success.
Question. Is anything known about $\Delta(s_\lambda)$ for $n>2$?