  The Virtual Element Method is well suited to the formulation of arbitrarily
regular Galerkin approximations of elliptic partial differential equations of
order $2p_1$, for any integer $p_1\geq 1$. In fact, the virtual element
paradigm provides a very effective design framework for conforming, finite
dimensional subspaces of $H^{p_2}(\Omega)$, $\Omega$ being the computational
domain and $p_2\geq p_1$ another suitable integer number. In this study, we
first present an abstract setting for such highly regular approximations and
discuss the mathematical details of how we can build conforming approximation
spaces with a global high-order continuity on $\Omega$. Then, we illustrate
specific examples in the case of second- and fourth-order partial differential
equations, that correspond to the cases $p_1=1$ and $2$, respectively. Finally,
we investigate numerically the effect on the approximation properties of the
conforming highly-regular method that results from different choices of the
degree of continuity of the underlying virtual element spaces and how different
stabilization strategies may impact on convergence.
