Conformal coupling of massive Klein-Gordon equation? Consider the massless Klein-Gordon equation of the form
$$\Big( \square_g - \xi R_g \Big) \phi \enspace = \enspace 0 \tag{1} \quad , $$
where $\square_g$ and $R_g$ are the d'Alembert operator and the Ricci-scalar-curvature with respect to the metric $g$ and $\xi$ is some constant. The case $\xi = 1/6$ is conformally invariant, which means that if $\phi$ solves equation (1), then under a transformation
$$ \tilde{g} \enspace \longmapsto \enspace \Omega^2 g \quad , \qquad \tilde{\phi} \enspace \longmapsto \enspace \Omega^{-1} \phi $$
the equation transforms as
$$ \Big( \square_g - \xi R_g \Big) \phi \enspace = \enspace \Omega^{-3} \Big( \square_{\tilde{g}} - \xi R_{\tilde{g}} \Big) \tilde{\phi} \enspace = \enspace 0 $$
and so the scalar-field $\tilde{\phi}$ solves the equation
$$ \Big( \square_{\tilde{g}} - \xi R_{\tilde{g}} \Big) \tilde{\phi} \enspace = \enspace 0.$$
Now the question is if there is a similar condition on $\xi$ (and/or $m^2$) such that the massive equation
$$\Big( \square_g - m^2 - \xi R_g \Big) \phi \enspace = \enspace 0 \tag{2} $$
becomes conformally invariant.

I suppose that since equation (1) with $\xi = 1/6$ is conformally invariant, that the additional terms emerging from a conformal transformation of the d'Alembert operator $\square_g$ cancel with the additional terms emerging from the conformal transformation of the $R_g$-term. I therefore furthermore suppose that adding this $m^2$-term destroys conformal invariance for all $m^2 \neq 0$. However, assuming $m^2$ e.g. not as a constant parameter but rather as a constant scalar-field that transforms just like the scalar-field $\phi$ would again reestablish conformal invariance. Is such an interpretation valid or is my original claim that there is no conformal invariance for $m^2 \neq 0$ true?
 A: Conformal invariance of this equation is not possible, if the mass $m$ is a constant scalar (conformal weight $0$). $m$ needs to have a nontrivial conformal weight and in fact, $\widetilde m=e^{\omega}m$ (conformal weight $1$) is the only possibility (even in dimensions other than $D=4$).
Theorem: Given the conformal transformations $\widetilde g_{\mu\nu}=e^{-2\omega}g_{\mu\nu}$ (conformal weight $-2$), $\widetilde m=e^{k\omega}m$ (conformal weight $k$) as well as $\widetilde\psi=e^{n\omega}\psi$ (conformal weight $n$), then the equation:
$$\left(\square-m^2-\xi R\right)\psi=0$$
is conformally invariant if and only if:
$$k=1,\qquad n=\frac{D-2}{2},\qquad \xi=\frac{D-2}{4(D-1)},$$
where $D$ is the dimension of spacetime.
So we have the remarkable result of no free variables despite the dimension $D$ as well as $k$ not even depending on it. (In the following proof, I have not included all the detailed calculations for clarity, but could provide them, if necessary.)
Proof: Consider the transformation of the upper equation under the given transformations, which will lead to conditions on $k$, $n$ and $\xi$ (showing the forward direction), which will simultaneously show, that those make the equation conformally invariant (showing the backwards direction).
The first term of the equation transforms as:
$$\widetilde\square\widetilde\psi
=e^{(n+2)\omega}\left(
\square\psi
+n\psi\square\omega
+(-D+2n+2)(\partial\omega)(\partial\psi)
+(-D+n+2)n\psi(\partial\omega)^2\right).$$
The term with $\partial\psi$ has to be eliminated putting $n=(D-2)/2$ as such a term cannot appear in the conformal transformations of the second or third term of the equation, as they don't contain derivatives of $\psi$, and hence can't be canceled to disappear otherwise. We get:
$$\widetilde\square\widetilde\psi
=e^{(n+2)\omega}\left(
\square\psi
+n(\square\omega
-n(\partial\omega)^2)\psi\right).$$
The third term of the equation transforms as:
$$\xi\widetilde{R}\widetilde\psi
=e^{(n+2)\omega}\xi\left(R+2(D-1)\left(\square\omega-\frac{(D-2)}{2}(\partial\omega)^2\right)\right)\psi.$$
This is a strong result as it confirmes $n=(D-2)/2$, which is the only possible choice that allows to cancel $\square\omega$ and $(\partial\omega)^2$ simultaneously with the equation before, which is only possible for:
$$\xi=\frac{n}{2(D-1)}
=\frac{D-2}{4(D-1)}.$$
We get:
$$\left(\widetilde\square-\xi\widetilde R\right)\widetilde\psi
=e^{(n+2)\omega}(\square-\xi R)\psi.$$
The second term of the equation transforms as:
$$\widetilde m^2\widetilde\psi
=e^{(n+2k)\omega}m^2\psi,$$
so $k=1$, when compared to the previous equation. $\square$
