  We devise and analyze two hybrid high-order (HHO) methods for the numerical
approximation of the biharmonic problem. The methods support polyhedral meshes,
rely on the primal formulation of the problem, and deliver $O(h^{k+1})$
$H^2$-error estimates when using polynomials of order $k\ge0$ to approximate
the normal derivative on the mesh (inter)faces. Both HHO methods hinge on a
stabilization in the spirit of Lehrenfeld--Sch\"oberl for second-order PDEs.
The cell unknowns are polynomials of order $(k+2)$ that can be eliminated
locally by means of static condensation. The face unknowns approximating the
trace of the solution on the mesh (inter)faces are polynomials of order $(k+1)$
in the first HHO method which is valid in dimension two and uses an original
stabilization involving the canonical hybrid finite element, and they are of
order $(k+2)$ for the second HHO method which is valid in arbitrary dimension
and uses only $L^2$-orthogonal projections in the stabilization. A comparative
discussion with the weak Galerkin methods from the literature is provided,
highlighting the close connections and the improvements proposed herein.
Additionally, we show how the two HHO methods can be combined with a
Nitsche-like boundary-penalty technique to weakly enforce the boundary
conditions. An originality in the devised Nitsche's technique is to avoid any
penalty parameter that must be large enough. Finally, numerical results
showcase the efficiency of the proposed methods and indicate that the HHO
methods can generally outperform discontinuous Galerkin methods and even be
competitive with $C^0$-interior penalty methods on triangular meshes.
