Let $a\_0$ and $b\_0$ be smooth compactly supported functions in $B \subset R^3$, $f\in C^1(\Omega)$, and define


$a\_n=f\Delta^{-1}(a\_{n-1})=-f(x)\int\_{B}a\_{n-1}(y)\Phi(x-y)dy$, $n\geq 1$


$b\_n=f\Delta^{-1}(b\_{n-1})=-f(x)\int\_{B}b\_{n-1}(y)\Phi(x-y)dy$, $n\geq 1$


where $\Phi$ is the fundamental solution of the laplacian. Suppose


$\int\_{B}(a\_n-b\_n)\varphi=0$ for every harmonic function $\varphi$ and every $n\geq 0$. Is $a\_0=b\_0$? When $f$ is a constant, the above means that $a\_0-b\_0$ is orthogonal to all functions $\varphi$ with $\Delta^n \varphi=0$ for some $n\in N$. This would imply $a\_0=b\_0$. Can this be generalized for a general coefficient $f$? 



Orthogonality to harmonic functions
This is an extended comment, not an answer. Suppose that $a\_0$ and $b\_0$ satisfy the conditions given in the statement of the problem and $u := b\_0-a\_0$ is not identically zero.


The function $u$ is orthogonal to harmonic functions. Since the difference between $\Phi(x-y)$ and the Green function $G\_\Omega(x,y)$ is a harmonic function, we have $\int\_\Omega \Phi(x-y) u(y) dy = \int\_\Omega G\_\Omega(x,y) u(y) dy =: G\_\Omega u(x)$. It follows that $f G\_\Omega u$, $f G\_\Omega(f G\_\Omega u)$ and so on are orthogonal to harmonic functions.


If $f$ is positive and bounded, the operator $u \mapsto G\_\Omega(f u)$, denoted by $G\_\Omega f$ for simplicity, is a compact, non-negative definite, self-adjoint operator on $L^2(\Omega, f(x)dx)$, and hence it has an orthonormal basis of eigenfunctions $g\_n$, with corresponding eigenvalues $\lambda\_n$ converging monotonically to zero.


Orthogonality of $(f G\_\Omega)^k u$ and a harmonic function $\phi$ is equivalent to orthogonality of $u$ and $(G\_\Omega f)^k \phi$. Write $f^{-1} u$ and $\phi$ in the basis $(g\_n)$: $f^{-1} u = \sum \alpha\_k g\_k$ and $\phi = \sum \beta\_k g\_k$. Then
$$
0 = \int\_\Omega u(x) (G\_\Omega f)^k \phi(x) = \sum \lambda\_n^k \alpha\_n \beta\_n .
$$
Let $M$ be the set of indices $m$ satisfying $\lambda\_m = \max\{\lambda\_n : \alpha\_n \ne 0\}$. Multiplying both sides of the above equality by $\lambda\_m^{-k}$ and passing to the limit as $k \to \infty$, we get
$$
0 = \sum\_{m \in M} \alpha\_n \beta\_n .
$$
Therefore, the function $v$ defined by $f^{-1} v = \sum\_{m \in M} \alpha\_n g\_n$ is orthogonal to all harmonic functions. Furthermore, $f^{-1} v$ is the eigenfunction of $G\_\Omega f$, that is, $G\_\Omega v(x) = \lambda\_m (f(x))^{-1} v(x)$.




---


This shows a different point of view on the original problem: take any function $v$ that is orthogonal to all harmonic functions, and evaluate $f(x) = v(x) / G\_\Omega v(x)$. Can $f$ be bounded?


This is not an equivalent formulation, as it does not tell anything about the support of $u$ (which was required to be a compact subset of $\Omega$). But it seems to be reasonably close; in particular, if we can show that $v$ with the above properties does not exist, then no solution $u$ to the original problem can exist.


I bet $f$ cannot be bounded, but I have no proof. A heuristic argument: for $x$ close to a boundary point $x\_0$ of $\Omega$, $G\_\Omega(x, y)$ is close to a multiple of the Poisson kernel $P\_\Omega(y, x\_0)$, which is harmonic, and therefore orthogonal to $v$. Therefore, $G\_\Omega v(x)$ *should* be of smaller order than $v(x)$ as $x$ converges to $x\_0$, and so $f$ cannot be bounded.




---


Edited: I forgot to include two relatively simple remarks.


(1) If $\Omega$ has no holes (that is, its complement is connected), then it does not matter whether we assume that $a\_n - b\_n$ is orthogonal to harmonic functions on $\Omega$ or on $\mathbb{R}^3$. Indeed, by [Runge's theorem for harmonic functions](http://onlinelibrary.wiley.com/doi/10.1112/jlms/s2-1.1.404/pdf), for any compact subset $K$ of $\Omega$, any harmonic function on $\Omega$ can be approximated uniformly on $K$ by harmonic polynomials.


(2) If $f$ is zero on some ball $B$, then the answer is **no**: any pair $a\_0, b\_0$ such that $a\_0 - b\_0$ is orthogonal to harmonic functions and $a\_0 - b\_0 = 0$ outside $B$ has the desired conditions. Indeed: $\int\_\Omega (a\_0(y) - b\_0(y)) \Phi(x - y) dy = 0$ outside $B$ because $\Phi(x - \cdot)$ is a harmonic function on (a neighbourhood of) $B$, and by remark (1), $a\_0 - b\_0$ is orthogonal to harmonic functions on (a neighbourhood of) $B$. This means that $a\_1 = b\_1$, and consequently $a\_n = b\_n$ for all $n \geqslant 1$. On the other hand, it is not difficult to construct a function $a\_0 - b\_0$ supported in $B$ and orthogonal to harmonic functions.


