Let $(M,g)$ be a closed (compact, no boundary) smooth $n$-dimensional Riemannian manifold. The Laplace–Beltrami operator $\Delta_g$ on $M$ has discrete spectrum $(\lambda_j)_j$ (indexed without multiplicity) with corresponding eigenfunctions $(\psi_j)_j$, normalized in $L^2(\rm{vol}_g)$. The MP zeta function of $\Delta_g$ at a point $p$ in $M$ is defined as $$\zeta^{\Delta_g}_p(s):=\sum_j \lambda_j^{-s} \psi_j(p)^2, \qquad s\in \mathbb C.$$

Question: Is something known about the continuity/differentiability of $p\mapsto \zeta^{\Delta_g}_p(s)$ as a real-valued function on $M$, say, for real $s>n/2$ ?

As far as I understand, the original paper by Minakshisundaram and Pleijel is only concerned with properties of $\zeta^{\Delta_g}_p(s)$ as a meromorphic function of $s$, holomorphic for $\Re(s)>n/2$. However, it seems to me that the authors do not discuss the properties of $p\mapsto \zeta^{\Delta_g}_p(s)$ for fixed $s$ (equivalently, of the polynomials $A$ in Eqn. (28) in their paper).

Based on Hörmander-type estimates on the $L^\infty$-norm of $\psi_j$ and Weyl's asymptotic, one can show that, for large enough $s$ (e.g. $s\geq n$), the function $p\mapsto \zeta^{\Delta_g}_p(s)$ is continuous (by showing the total convergence of the series). This is however unsatisfactory, for I would hope that $p\mapsto \zeta^{\Delta_g}_p(s)$ is continuous (possibly even smooth) for all $s>n/2$, which is the threshold originating in Weyl's asymptotic for the $\lambda_j$'s. Is there any result in this direction which holds for all $s>n/2$?