Let $(M,g)$ be a closed (compact and without boundary) and oriented Riemannian manifold and let us consider the Poisson equation for a smooth function $\varphi$:

$\Delta \phi = f$,

where $f$ is a fixed smooth function on $M$. It is then well-known that the previous Poisson equation has a solution if and only if

$\int_M f = 0$,

and in that case, the solution is unique up to the addition of a constant. I was wondering if this result generalizes to a system of Poisson equations for smooth functions $\phi_1$ and $\phi_2$:

$\Delta \phi_1 = f_1 \phi_1 +f_2 \phi_2\, , \qquad \Delta \phi_2 = g_1 \phi_1 + g_2 \phi_2$

for fixed functions $f_1$, $g_1$, $f_2$ and $g_2$. Clearly, a necessary condition for $\phi_1$ and $\phi_2$ to be solutions is:

$ \int_M (f_1 \phi_1 +f_2 \phi_2) = 0\, , \qquad \int_M (g_1 \phi_1 + g_2 \phi_2) = 0$

I wonder if the converse is also true. I am sure this problem has been extensively studied, but I have not found a reference about it.

Thanks.