An $hp$ version of interface penalty finite element method ($hp$-IPFEM) is
proposed for elliptic interface problems in two and three dimensions on
unfitted meshes. Error estimates in broken $H^1$ norm, which are optimal with
respect to $h$ and suboptimal with respect to $p$ by half an order of $p$, are
derived. Both symmetric and non-symmetric IPFEM are considered. Error estimates
in $L^2$ norm are proved by the duality argument.