  In this paper, we prove that nonnegative polyharmonic functions on the upper
half space satisfying a conformally invariant nonlinear boundary condition have
to be the "\emph{polynomials} plus \emph{bubbles}" form. The nonlinear problem
is motivated by the recent studies of boundary GJMS operators and the
$Q$-curvature in conformal geometry. The result implies that in the conformal
class of the unit Euclidean ball there exist metrics with a single singular
boundary point which have flat $Q$-curvature and constant boundary
$Q$-curvature. Moreover, all of such metrics are classified. This phenomenon
differs from that of boundary singular metrics which have flat scalar curvature
and constant mean curvature, where the singular set contains at least two
points. A crucial ingredient of the proof is developing an approach to separate
the higher order linear effect and the boundary nonlinear effect so that the
kernels of the nonlinear problem are captured.
