Modified mean value property Let $L=\Delta + c$ in 3 dimensions, where $c$ is a positive constant.
I met this modified mean value property of a solution $u$ of $Lu=0$ as 
$$u(\xi)=\frac{\sqrt{c}\rho}{sin(\sqrt{c}\rho)}\frac{1}{4\pi \rho^2}\int_{\partial B(\xi,\rho)} u(x)d\sigma(x)$$
where $sin(\sqrt{c}\rho) \ne 0$, and $\sigma$ denotes the surface measure.
I'm trying to prove it. The following is what I did.
First I figured out that $$K(x,\xi)=-\frac{cos(\sqrt{c}r)}{4\pi r}, r=|x-\xi|$$ is a fundamental solution of $L$ with pole $\xi$. Then I tried to apply the second Green identity (which is also true for $L$ for sure)
$$\int_{\Omega}(LG) u-\int_{\Omega}G (Lu)=\int_{\partial \Omega} \frac{\partial G}{\partial \nu} u- G \frac{\partial u}{\partial \nu}$$ For fixed $\rho$, let $G=K+\frac{cos(\sqrt{c}\rho)}{4\pi \rho}$ and $\Omega=B_{\rho}(\xi)$, so that $G$ vanishes on the boundary of $\Omega$, and $LG=\delta_{\xi}+c\frac{cos(\sqrt{c}\rho)}{4\pi \rho}$, where $\delta_{\xi}$ denotes the Dirac mass at $\xi$. Then the identity above becomes $$u(\xi)+c\frac{cos(\sqrt{c}\rho)}{4\pi \rho}\int_{B_{\rho}(\xi)}u(x)dx=\frac{\sqrt{c}\rho sin(\sqrt{c}\rho)+cos(\sqrt{c}\rho)}{4 \pi \rho^2} \int_{\partial B(\xi,\rho)} u(x)d\sigma(x)$$
Then no matter how hard I tried, I just could not get the desired form. 
I also tried to take the derivative with respect to $\rho$ of the desired form, and I got a complicated form which doesn't seem to be zero. Then I took the second derivative, and things got worse...
Can anyone either give the right reference about the solution of this operator $L$, or help me figure out how to obtain the desired form? Thanks!
 A: Ah, it is one of these cases when you better know your formula in all dimensions. 
    It is 
    $$
  \int_{|x-y|=\rho} u(x) d\sigma(x) = (2\pi)^{n/2} \cdot \frac{J_{n/2-1}(\sqrt{c}\rho)}{(\sqrt{c}\rho)^{n/2-1}} \cdot u(y)
 $$
    One of the ways to prove it:
    $$
 \Delta u + c u = 0,
 $$
    $$
 \int_{|x-y|\leq R} \Delta u(x) + c u(x)  dx = 0.
 $$
    Apply Green's formula to the first integrand, 
    $$
 \int_{|x-y|= R} \frac{\partial u}{\partial r}(x) d\sigma(x) + \int_{|x-y|\leq R} c u(x)  dx = 0,
 $$
    (here $\frac{\partial u}{\partial r}(x) = \nabla u(x)\cdot x/|x|$).
    $$
 R^{n-1}\int_{|x|=1} \frac{\partial u}{\partial r}(y-Rx) d\sigma + c \int_{|x-y|\leq R} u(x) dx = 0.
 $$
    Differentiate w.r. to $R$ 
    $$
 R^{n-1}\int_{|x|=1} \frac{\partial^2 u}{\partial r^2}(y-Rx) d\sigma +(n-1)R^{n-2} \int_{|x|=1} \frac{\partial u}{\partial r}(y-Rx) d\sigma + c R^{n-1} \int_{|x|=1} u(y-Rx) d\sigma = 0,
 $$
    $$
 \int_{|x|=1} \frac{\partial^2 u}{\partial r^2}(y-Rx) d\sigma +\frac{n-1}{R}\int_{|x|=1} \frac{\partial u}{\partial r}(y-Rx) d\sigma + c \int_{|x|=1} u(y-Rx) d\sigma = 0.
 $$
    $$
 \frac{d^2}{dR^2}\int_{|x|=1} u (y-Rx) d\sigma +\frac{n-1}{R} \frac{d}{dR}\int_{|x|=1} u(y-Rx) d\sigma + c \int_{|x|=1} u(y-Rx) d\sigma = 0.
 $$
    So $v(R) = \int_{|x|=1} u (y-Rx) d\sigma $ is a solution of 
    $$
 v''(R) + \frac{n-1}{R} v'(R) + c v(R)=0.
 $$
    This equation is one of the forms of Bessel's equation 
    $$
 x^2 y''+ xy' + (x^2-\alpha^2) y = 0. 
 $$
    In particular, it is not difficult to show that 
    $$
 v(R) = C_1 \frac{J_{n/2-1}(R\sqrt{c})}{R^{n/2-1}},
 $$
    here $J_\nu$ is the Bessel function of the first kind.
    To get the right constant just send $R\to 0$: on the one hand 
    $$
 \lim_{R\to 0} v(R)  = \omega_{n-1} \cdot u(y),
 $$ 
    on the other 
    $$
 \lim_{R\to 0} v(R) = C_1\lim_{R\to 0} \frac{J_{n/2-1}(R\sqrt{c})}{R^{n/2-1}} = C_1 \frac{(\sqrt{c}/2)^{n/2-1}}{2\pi^{n/2}} \omega_{n-1}. 
 $$
    So $C_1 = (2\pi)^{n/2} c^{1/2-n/4} u(y)$.
