Let $\Omega$ be a bounded domain in $\mathbb{R}^3$ and $f\_1,f\_2 \in C^2(\bar{\Omega})$. Suppose 


$\int\_{\Omega}(f\_2-f\_1)\varphi \, dx=0$ and $\int\_{\Omega}(f\_2 \Delta^{-1} f\_2- f\_1 \Delta^{-1} f\_1)\varphi \, dx =0$


for all harmonic functions $\varphi$, where $\Delta^{-1}f=\int\_{\Omega}f(y)\Phi(x-y) \, dy$ and $\Phi$ is the fundamental solution of the Laplacian. Is $f\_1=f\_2$? 



Functions orthogonal to harmonic functions
Assume without loss of generality that $0\in \Omega$. 


Let $f\_1, f\_2$ be radial, with support contained in $\Omega$. 


By the mean value property of harmonic functions, you have


$$ \int f\_2 \varphi ~\mathrm{d}x = 4\pi \varphi(0) \int\_0^{\infty} r^2 f\_2(r) ~\mathrm{d}r = \varphi(0) \int f\_2(r) ~\mathrm{d}x. $$ 


Now let $f\_1, f\_2$ be arbitrary non-trivial radial functions such that $\int f\_1 ~\mathrm{d}x = \int f\_2 ~\mathrm{d}x= 0$.


As $\Delta^{-1} f\_1$ and $\Delta^{-1}f\_2$ are also radial, we have that $$\int f\_1 \Delta^{-1} f\_1 \varphi ~\mathrm{d}x = \varphi(0) \int f\_1 \Delta^{-1} f\_1 ~\mathrm{d}x$$
whenever $\varphi$ is harmonic. So then by replacing $f\_1 \mapsto \lambda\_1 f\_1$ and $f\_2 \mapsto \lambda\_2 f\_2$ we can certainly arrange for 


$$ \int (\lambda\_1)^2 f\_1 \Delta^{-1} f\_1 - (\lambda\_2)^2 f\_2 \Delta^{-2} f\_2 ~\mathrm{d}x = 0 $$


and so as long as our seeds $f\_1$ and $f\_2$ have, for example, unequal supports, we arrive at a counterexample. 




---


In the comments below the OP asked for examples where $f\_1$ and $f\_2$ are signed. 
Hopefully the computations below are (more or less) correct:


Roughly speaking, let $1 > r > \rho$ be fixed numbers. The idea is to let $f\_1$ be the sum of the surface measure of the sphere of radius 1 with $a\_\rho$ times the surface measure of the sphere of radius $\rho$, while $f\_2$ is $a\_r$ times the surface measure of the sphere of radius $r$. Then the first condition requires


$$ 1 + a\_\rho \rho^2 = a\_r r^2 $$


For the surface measure of a sphere of radius $\rho$, $\Delta^{-1}$ can be computed to be the function equaling $C |x|^{-1} \rho^2$ for some universal constant $C$. So you would have that the condition on $f\_1 \Delta^{-1} f\_1$ to give rise to the condition


$$ 1 + a\_\rho^2 \rho + 2 a\_\rho \rho^2 = a\_r^2 r $$


$$ \implies 1 + a\_\rho^2 \rho + 2 a\_\rho \rho^2 = r^{-3} (1 + 2 a\_\rho \rho^2 + a\_\rho^2 \rho^4) $$


$$ \implies 1 - r^3 + 2 a\_\rho \rho^2 (1 - r^3) + a\_\rho^2 \rho (\rho^3 - r^3) = 0$$


Since $1 > r > \rho$ this last equation has a positive root $a\_\rho$, which then gives a positive solution $a\_r$. 


To move from measures to functions: instead of using surface measures use $C^\infty\_c$ functions highly concentrated near the corresponding surfaces. The computations would be slightly less explicit but the relationships between the coefficients $a\_r, a\_\rho$ would still hold up to a (small) perturbation, and you still can find a counterexample. 


