Particle-In-Cell codes are widely used for plasma physics simulations. It is often the case that particles within a computational cell need to be split to improve the statistics or, in the case of non-uniform meshes, to avoid the development of fictitious self-forces. Existing particle splitting methods are largely empirical and their accuracy in preserving the distribution function has not been evaluated in a quantitative way. Here we present a new method specifically designed for codes using adaptive mesh refinement. Although we point out that an exact, distribution function preserving method does exist, it requires a large number of split particles and its practical use is limited. We derive instead a method that minimizes the cost function representing the distance between the assignment function of the original particle and that of the sum of split particles. Depending on the interpolation degree and on the dimension of the problem, we provide tabulated results for the weight and position of the split particles. This strategy represents no overhead in computing time and for a large enough number of split-particles it asymptotically tends to the exact solution.