This paper presents a grid-free simulation algorithm for the fully three-dimensional Vlasov--Poisson system for collisionless electron plasmas. We employ a standard particle method for the numerical approximation of the distribution function. Whereas the advection of the particles is grid-free by its very nature, the computation of the acceleration involves the solution of the nonlocal Poisson equation. To circumvent a volume mesh, we utilize the fast boundary element method, which reduces the three-dimensional Poisson equation to a system of linear equations on its two-dimensional boundary. This gives rise to fully populated matrices which are approximated by the $\mathcal H^2$-technique, reducing the computational time from quadratic to linear complexity. The approximation scheme based on interpolation has been shown to be robust and flexible, allowing a straightforward generalization to vector-valued functions. In particular, the Coulomb forces acting on the particles are computed in linear complexity. In first numerical tests, we validate our approach with the help of classical nonlinear plasma phenomena. Furthermore, we show that our method is able to simulate electron plasmas in complex three-dimensional domains with mixed boundary conditions in linear complexity.