  The realization of a standard Adaptive Finite Element Method (AFEM) preserves
the mesh conformity by performing a completion step in the refinement loop: in
addition to elements marked for refinement due to their contribution to the
global error estimator, other elements are refined. In the new perspective
opened by the introduction of Virtual Element Methods (VEM), elements with
hanging nodes can be viewed as polygons with aligned edges, carrying virtual
functions together with standard polynomial functions. The potential advantage
is that all activated degrees of freedom are motivated by error reduction, not
just by geometric reasons. This point of view is at the basis of the paper [L.
Beirao da Veiga et al., Adaptive VEM: stabilization-free a posteriori error
analysis and contraction property, SIAM Journal on Numerical Analysis, vol. 61,
2023], devoted to the convergence analysis of an adaptive VEM generated by the
successive newest-vertex bisections of triangular elements without applying
completion, in the lowest-order case (polynomial degree k=1).
  The purpose of this paper is to extend these results to the case of VEMs of
order k>1 built on triangular meshes. The problem at hand is a
variable-coefficient, second-order self-adjoint elliptic equation with
Dirichlet boundary conditions; the data of the problem are assumed to be
piecewise polynomials of degree k-1. By extending the concept of global index
of a hanging node, under an admissibility assumption of the mesh, we derive a
stabilization-free a posteriori error estimator. This is the sum of
residual-type terms and certain virtual inconsistency terms (which vanish for
k=1). We define an adaptive VEM of order k based on this estimator, and we
prove its convergence by establishing a contraction result for a linear
combination of (squared) energy norm of the error, residual estimator, and
virtual inconsistency estimator.
