  We construct and analyze conformally invariant random fields on 4-dimensional
Riemannian manifolds $(M,g)$. These centered Gaussian fields $h$, called
\emph{co-biharmonic Gaussian fields}, are characterized by their covariance
kernels $k$ defined as the integral kernel for the inverse of the \emph{Paneitz
operator} \begin{equation*}\mathsf p=\frac1{8\pi^2}\bigg[\Delta^2+
  \mathsf{div}\left(2\mathsf{Ric}-\frac23\mathsf{scal}\right)\nabla \bigg].
\end{equation*} The kernel $k$ is invariant (modulo additive corrections) under
conformal transformations, and it exhibits a precise logarithmic divergence
$$\Big|k(x,y)-\log\frac1{d(x,y)}\Big|\le C.$$ In terms of the co-biharmonic
Gaussian field $h$, we define the \emph{quantum Liouville measure}, a random
measure on $M$, heuristically given as \begin{equation*}
  d\mu(x):= e^{\gamma h(x)-\frac{\gamma^2}2k(x,x)}\,d \text{vol}_g(x)\,,
\end{equation*} and rigorously obtained a.s.~for $|\gamma|<\sqrt8$ as weak
limit of the RHS with $h$ replaced by suitable regular approximations
$(h_\ell)_{\ell\in\mathbb N}$.
  For the flat torus $M=\mathbb T^4$, we provide discrete approximations of the
Gaussian field and of the Liouville measures in terms of semi-discrete random
objects, based on Gaussian random variables on the discrete torus and piecewise
constant functions in the isotropic Haar system.
