Solution (source code)

= Solution

In the paper's convention, $\dot h_{ij}=2NK_{ij}+D_iN_j+D_jN_i$. The <Legendre transform in mechanics> therefore gives
$$
H=\int d^3x\left[\sqrt hN(K_{ij}K^{ij}-K^2-{}^{(3)}R)+2\pi^{ij}D_iN_j\right].
$$
Indeed $\pi^{ij}K_{ij}=\sqrt h(K_{ij}K^{ij}-K^2)$, and inversion of the <canonical momentum of the spatial metric> gives
$$
K_{ij}K^{ij}-K^2=\frac1h\left(\pi^{ij}\pi_{ij}-\frac12\pi^2\right).
$$
To integrate the shift term without ambiguity about <tensor densities>, put $p^{ij}=\pi^{ij}/\sqrt h$. Metric compatibility of the <spatial covariant derivative> gives, up to a boundary term,
$$
\int d^3x\,2\sqrt h\,p^{ij}D_iN_j
=-\int d^3x\,2\sqrt h\,N^iD_jp^j{}_i.
$$
Thus the <Hamiltonian formulation of general relativity> has
$$
\boxed{\begin{aligned}
H&=\int d^3x\sqrt h(N\mathcal H+N^i\mathcal H_i),\\
\mathcal H&=\frac{\pi^{ij}\pi_{ij}-\pi^2/2}{h}-{}^{(3)}R,\\
\mathcal H_i&=-2D_j\left(\frac{\pi^j{}_i}{\sqrt h}\right).
\end{aligned}}
$$
Variation of the <lapse function> imposes the <Hamiltonian constraint> $\mathcal H=0$, while variation of the <shift vector> imposes the <momentum constraint> $\mathcal H_i=0$.