I am working on Q-curvature:
$$ Q^n_g= \frac{1}{2(n-1)} \Delta_g R_g+\frac{n^3-4n^2+16n-16}{8(n-1)^2(n-2)^2}\:R_g-\frac{2}{(n-2)^2}\:| Ric |^2_g $$
See, for example, the expression (3) in https://arxiv.org/pdf/0803.4331.pdf. However, the references I encounter only presents the expression of Q-curvature, but they do not indicate how to get it. I would like to know how to get it. Thanks in advance for some collaboration.
The motivation for this formula is that the conformal linearization of the $Q$-curvature is given by the Paneitz operator, a specific formally self-adjoint, conformally covariant operator with leading-order term a power of the Laplacian. So your question should really be "how do I get the formula for the Paneitz operator". Branson gave a simple derivation of this as Theorem 1.21 in his article "Differential operators canonically associated to a conformal structure". For convenience, I will sketch the method.
Let $$ P = \frac{1}{n-2} \left( \mathrm{Ric} - Jg \right) $$ denote the Schouten tensor, where $J=\mathrm{tr}_g P$ is the trace of the Schouten tensor. Your formula for $Q$ is equivalently written $$ Q = -\Delta J - 2\lvert P\rvert^2 + \frac{n}{2}J^2 . $$ The conformal transformation law for the Ricci tensor implies that $$ \left. \frac{\partial}{\partial t}\right|_{t=0} P^{e^{2tu}g} = -\nabla^2u , $$ where $\nabla^2u$ denotes the Hessian of $u$. The conformal transformation law for the Laplacian implies that $$ \left. \frac{\partial}{\partial t}\right|_{t=0} e^{t(2-w)u} \Delta_{e^{2tu}g} (e^{twu}f) = (n+2w-2)\langle\nabla u,\nabla f\rangle + wf\Delta u $$ for any $w\in\mathbb{R}$. Putting these together (and a similar formula for the conformal variation of the divergence $\delta$) implies that the Paneitz operator $$ Lf := \Delta^2 f + \delta \left( (4P - (n-2)Jg)(\nabla f) \right) + \frac{n-4}{4}Qf $$ is infinitesimally conformally covariant; i.e. $$ \left. \frac{\partial}{\partial t}\right|_{t=0} e^{\frac{n+4}{2}tu} L_{e^{2tu}g} \left( e^{-\frac{n-4}{2}tu} f \right) = 0 $$ for all smooth functions $u$ and $f$ and all metrics $g$. This implies that the Paneitz operator is conformally invariant; i.e. $$ e^{\frac{n+4}{2}u} L_{e^{2u}g} \left( e^{-\frac{n-4}{2}u}f \right) = L_g(f) $$ for all smooth functions $u$ and $f$.