  In this work, we propose a novel formulation for the solution of partial
differential equations using finite element methods on unfitted meshes. The
proposed formulation relies on the discrete extension operator proposed in the
aggregated finite element method. This formulation is robust with respect to
the location of the boundary/interface within the cell. One can prove enhanced
stability results, not only on the physical domain, but on the whole active
mesh. However, the stability constants grow exponentially with the polynomial
order being used, since the underlying extension operators are defined via
extrapolation. To address this issue, we introduce a new variant of aggregated
finite elements, in which the extension in the physical domain is an
interpolation for polynomials of order higher than two. As a result, the
stability constants only grow at a polynomial rate with the order of
approximation. We demonstrate that this approach enables robust high-order
approximations with the aggregated finite element method. The proposed method
is consistent, optimally convergent, and with a condition number that scales
optimally for high order approximation.
