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.

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, 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.