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$).