I am looking to numerically calulate the diagonal of Green's function. I am interested in Green's functions of elliptic PDEs and in those that arise from stochastic processes (discrete and continuous). I am aware of the connection - that is why I am interested in whatever I can find. Most literature I found is concerned with showing the diagonal is not regular in some sense and with bounds on it.

Let $G$ be the Green's function of $(-\Delta + \kappa)^2$ ($\kappa$ a constant killing rate) on $\Omega \subseteq \mathbb{R}^d, d \leq 3$. This operator needs two boundary conditions on $\Omega$. Take these conditions to be any of the standard three homogeneous boundary conditions taken with each application of the operator.

Is there a PDE that $G(x,x)$ satisfies? The probabilistic interpretation of the of Green's function is the expected number of visits to the start point (at least for random walks on graphs). I expect the diagonal to be harmonic if $k=0$. Is that true? Reasonable?

Here's what I've tried: For $d=1$, taking $y(x) = x$ and differentiating $\frac{d^2 G(x,y(x))}{dx^2}$, I get a term $\frac{\partial^2 G}{\partial x \partial y}|_{y=x}$ which I don't know what to do with. Using Mercer's theorem (assuming eigenfunctions are real, at least for now):

$$ G(x,y) = \sum_n \lambda_{n} \Psi_{n}(x) \Psi_n(y) $$

I get

$$ \frac{d^2 G(x,x)}{dx^2} = \sum_n \lambda_{n} \frac{d^2 \Psi_n^2(x)}{dx^2}\\ % % =2\sum_n \lambda_{n} \bigg (\Psi_n''(x)\Psi_n(x) + [\Psi_n'(x)]^2 \bigg )\\ % % =2\frac{\partial ^2 G(x,y(x))}{\partial x^2} + \sum_n \lambda_{n} [\Psi_n'(x)]^2\\ $$

which is (not surprisingly) pretty much the same problem. I considered taking an eigenvalue expansion of $(-\Delta + \kappa)^{-1}$ but I'm hoping for a better approach.

References I know: The authors of Diagonals of Green's functions and applications consider Green's functions for Sturm-Liouville operators that look similar to the operator I described above. They consider only one dimension - time. I don't know if this can be genralized to more than one spatial dimension. They do present a differential equation for the diagonal of the corresponding Green's function. However, their operator has non-homogeneous Dirichlet BC. They say their derivation was originally proposed in Semilinear Elliptic Equations Involving Critical Sobolev Exponents but I couldn't find any reference to diagonals in the latter.