Double-layer potentials on Riemannian manifolds Let $M$ be a compact Riemannian manifold, and let $S \subset M$ be a smooth hypersurface which divides $M$ into two domains $D_1$, $D_2$. Let also $g \colon S \to \mathbb R$ be a smooth function (everything is, say, $C^\infty$).
Is it true that there always exist (unique up to a common additive constant) harmonic functions $U_1 \colon D_1 \to \mathbb R$ and $U_2 \colon D_2 \to \mathbb R$ which are smooth up to $S$ and such that on $S$ we have:
1) $U_1 - U_2 = g$,
2) $\nabla U_1 - \nabla U_2$ is tangent to $S$ (i.e. the normal derivatives of $U_1$ and $U_2$ on $S$ coincide)?
This problem may be reformulated as a Poisson equation $\Delta U = \phi$ for a single function $U$, where $\phi$ is a distribution defined by
$
\langle \phi , f \rangle = \int_S g \cdot i_{\nabla f}\mu\,. 
$
Here we assume that $S$ is oriented in such a way that $\partial D_1 \ = S$, $\partial D_2 = -S$, and $\mu$ is the Riemannian volume form.
Then one can try to solve the problem by computing the convolution of $\phi$ with Green's function. The only problem will be to show that the resulting expression will have the regularity that we need. I believe this should not be too difficult (if true), but the whole problem sounds so natural to me that I am sure this must be well known.
UPD: When $M$ is a domain in $\mathbb R^n$, the solution to the above problem is given by the double-layer potential. The latter is the convolution of $\phi$ defined above with Green's function. However, I am still looking for a reference explaining this for a Riemannian manifold (in this case, the double-layer potential can be defined by the same formula as in the flat case; the only problem is to show that it has the required behavior at $S$; I expect that this should be so, because the asymptotic behavior of Green's function near the singularity is the same as in $\mathbb R^n$).
 A: Here is another approach : ([cf. Taylor PDE II Qualitative Studies of Linear Equations, chapter 7, section 11] 
We introduce the Dirichlet to Neumann operators 
 $N_i\colon \mathcal{C}^\infty(S)\rightarrow \mathcal{C}^\infty(S)$ defined by
 $$N_i(g)=\left.\frac{\partial}{\partial \nu_i}\right|_{\Sigma}h_i$$
 where $h_i$ is the harmonic extension of $g$ on $D_i$ and $\vec{\nu_i}$ is the unit outward normal to $S\subset D_i$. And $N=N_1+N_2$.
 These $3$ operators are self-adjoint Pseudo-differential elliptic operators of order $1$ and 
 $$\ker N_1=\ker N_2=\ker N=\{Constants\}.$$
You need to find two functions $g_1,g_2 \colon S \rightarrow \mathbb{R}$ such that
 $$  \left\{\begin{array}{l} g_1-g_2=g \\ N_1 (g_1)+N_2(g_2)=0 \end{array}\right.$$
Hence you need to find
$g_2$ such that 
$$N(g_2)=-N_1(g).$$
The fact that $N_1$ and $N$ are self-adjoint Pseudo-differential elliptic operators of order $1$ with the same kernel tells us that
this equation is solvable and that if $g\in H^{s}(S)$ then $g_1,g_2\in H^{s}(S)$ and the harmonic extensions of $g_1$ and $g_2$ are $H^{s+\frac12}(D_i)$.
A: At least as a sort of superficial/place-holder response... in the compact Riemannian manifold situation, the Laplace-Beltrami operator has a very nice resolvent. In particular, the corresponding "abstract" global Sobolev spaces have the property that the transition maps are at least Hilbert-Schmidt, and so on. A not-so-subtle, but global (i.e., not purely local...) consideration shows that any differential equation such as you pose can be solved by eigenfunction expansions. (Pointwise convergence of such expansions is a relatively subtle thing, but $L^2$ and global Sobolev convergence is/are "gimme"s.) 
