Showing that the conformal Laplacian is a conformally covariant operator 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!
 A: [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
