Consider a Riemannian manifold of dimension $n\geq3$. Consider the conformal Laplacian \begin{equation*} P_g=\Delta_g+\frac{n-2}{4(n-1)}R_g, \end{equation*} where $R_g$ is the scalar curvature associated to the metric $g$.
I would like to show that $P_g$ is a conformally covariant operator in the sense that under the conformal change $\tilde{g}=e^{2f}g$, the following transformation law is satisfied: \begin{align*} P_{\tilde{g}}&=e^{-(\frac{n}{2}+1)f}P_g \,e^{(\frac{n}{2}-1)f}\\ &= e^{-(\frac{n}{2}+1)f}\Delta_g \,e^{(\frac{n}{2}-1)f}+\frac{n-2}{4(n-1)}R_g\,e^{-2f} \tag{1} \end{align*} I know that under such a conformal change, the Laplace-Beltrami operator transforms like \begin{equation*} \Delta_{\tilde{g}}=e^{-2f}\Delta_g-(n-2)e^{-2f}g^{ij}\frac{\partial f}{\partial x_j}\frac{\partial}{\partial x_i}. \end{equation*} For the scalar curvature, one can make the substitution $e^{2f}=\varphi^{4/(n-2)}$ (where $\varphi$ is positive) to get \begin{equation*} R_{\tilde{g}}=\varphi^{-(n+2)/(n-2)}\bigg(4\frac{n-1}{n-2}\Delta_g \varphi + R_g\varphi\bigg) \end{equation*} which is really just \begin{equation*} R_{\tilde{g}}=4\frac{n-1}{n-2}e^{-(\frac{n}{2}+1)f}\Delta_g \, e^{(\frac{n}{2}-1)f}+R_g e^{-2f}. \end{equation*} So, I get \begin{equation*} \tag{2} P_{\tilde{g}}= e^{-2f}\Delta_g-(n-2)e^{-2f}g^{ij}\frac{\partial f}{\partial x_j}\frac{\partial}{\partial x_i} + e^{-(\frac{n}{2}+1)f}\Delta_g \, e^{(\frac{n}{2}-1)f}+\frac{n-2}{4(n-1)}R_g\, e^{-2f}. \end{equation*}
The problem is that I do not see how (1) is the same as (2). Have I made a mistake somewhere? Or do the first two terms in (2) somehow cancel each other out?
Any help would be greatly appreciated!
[An incomplete draft]
I would adjust the notation to make it more suitable for the upcoming calculation. The conformally rescaled metric will be denoted by $\widehat{g} = e^{2\varphi} g$, where $g$ is a Riemannian metric on a given manifold $M$ (the consideration is purely local, so we can be a little less specific here). It is customary then to denote the quantities corresponding to the conformally rescaled metric $\widehat{g}$ by placing the $\widehat{}$ on top of the respective symbol. Thus, the Levi-Civita connection $\nabla^{\widehat{g}}$ of the rescaled metric $\widehat{g}$ will be denoted simply by $\widehat{\nabla}$, and $\nabla$ will mean the Levi-Civita connection, corresponding to the metric $g$. Furthermore, $\Delta_{\widehat{g}}$ will be denoted by $\widehat{\Delta}$, whereas $\Delta_{g}$ will be simply denoted as $\Delta$. The same conventions apply to the scalar curvatures $\widehat{R}$ and $R$, and so on.
The reason why I make all these preparation is that we need to recall a few formulas for the conformal rescaling of the Levi-Civita connection. We are going to use this formula for the case of $1$-forms: $$ \widehat{\nabla}_i \omega_j = \nabla_i \omega_j - \omega_i \varphi_j - \varphi_i \omega_j + \varphi^k \omega_k g_{i j} $$ where $\varphi_i := \nabla_i \varphi$.
Using this formula, it is straightforward to show that $$ \widehat{\Delta} f = e^{-2 \varphi} \big( \Delta f - (n - 2) \varphi^k \nabla_k f \big) $$ and, with a little bit more work, that the scalar curvature rescales as (see. e.g. here) $$ \widehat{R} = e^{-2 \varphi} \big( R + 2(n - 1) \Delta \varphi - (n - 2)(n - 1) \varphi^k \varphi_k \big) $$ It is not too difficult to convince yourself that there is no linear combination of operators $\Delta$ and $R$ (regarded as acting by multiplication) such that all the extra terms in their transformation would cancel out.
The trick is to consider the action of $\Delta$ and $R$ on weighted functions, that is functions of the form $e^{a \varphi} f$.
One can verify that $$ \widehat{\nabla}_i (e^{a \varphi} f) = e^{a \varphi} \big( \nabla_i f + a \varphi_i f \big) $$ and, using this, that $$ \widehat{\Delta} (e^{a \varphi} f) = e^{(a - 2) \varphi} \bigg( \Delta f + (n + 2 a -2) \varphi^k \nabla_k f + a \big( \nabla^k \varphi_k + (n + a -2) \varphi^k \varphi_k \big) f \bigg) $$
For $a = 1 - n/2$ we get $$ \widehat{\Delta} (e^{ (1 - n/2) \varphi} f) = e^{(- 1 - n/2) \varphi} \bigg( \Delta f - \frac{n-2}{4} \big( 2 \Delta \varphi + (n - 2) \varphi^k \varphi_k f \bigg) $$
At the same time $$ \widehat{R}(e^{ (1 - n/2) \varphi} f) = e^{(- 1 - n/2) \varphi} \bigg( R - (n - 1) \big( 2 \Delta \varphi + (n - 2) \varphi^k \varphi_k \big) f \bigg) $$
At this point it should be easy to see that the operator $P_g$, as defined in the question, exhibits the necessary covariance property.
TODO: check the signs.
References.
Jan Slovak, Natural Operators on Conformal Manifolds, http://www.math.muni.cz/~slovak/Papers/habil.pdf
M.Eastwood, Conformally invariant differential operators, https://maths-people.anu.edu.au/~eastwood/fayetteville2.pdf
List of formulas in Riemannian geometry, https://en.wikipedia.org/wiki/List_of_formulas_in_Riemannian_geometry