The detection and characterization of debris in the integral-of-motion space is a promising avenue to uncover the hierarchical formation of the Milky Way. Yet, the fact that the integrals do not remain constant during the assembly process adds considerable complexity to this approach. Indeed, in time-dependent potentials tidal substructures tend to be effaced from the integral-of-motion space through an orbital diffusion process, which naturally leads to the formation of a ‘smooth’ stellar halo. In this talk I will introduce a new probability theory that describes the evolution of collisionless systems subject to a time-dependent potential. The new theory can be used to reconstruct the hierarchical assembly of our Galaxy through modelling the observed distribution of accreted stars in the integral-of-motion space.