[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
The issue is that the first formula comes from taking components of $X$ with respect to the orthonormal basis $e_r,e_\theta$, whereas the differential geometry formula is in terms of components with respect to the coordinate basis $\partial/\partial r, \partial/\partial\theta$. In particular, the unit vector $e_\theta = \frac1r\partial/\partial\theta$. So the $X^\theta$ appearing in the two formulas are different.
Best Answer
Here's an relatively straightfoward, coordinate-free way to approach the problem:
The divergence of a metric $g$ actually depends only on its volume form $\mu$ (for either orientation, it doesn't matter which). It can be characterized by $$\mathcal{L}_X \mu = (\text{div } X) \mu.$$
Now, the volume form $\widehat{\mu}$ of a conformally related metric $\widehat{g} = e^{2f} g$ is $$\widehat{\mu} = e^{n f} \mu;$$ substituting in the above characterization of the divergence---and using the Leibniz identity $$\mathcal L_X (h \alpha) = (X \cdot h) \alpha + h \mathcal L_X \alpha$$ for a function $h$ and form $\alpha$---gives an identity relating the divergences of the two metrics.