A theoretical framework is developed to describe, in the limit of small but finite Re, the evolution of dilute clusters of sedimenting particles. Here, Re =aU/ν is the particle Reynolds number, where a is the radius of the spherical particle, U its settling velocity, and ν the kinematic viscosity of the suspending fluid. The theory assumes the disturbance velocity field at sufficiently large distances from a sedimenting particle, even at small Re, to possess the familiar source--sink character; that is, the momentum defect brought in via a narrow wake behind the particle is convected radially outwards in the remaining directions. It is then argued that for spherical clusters with sufficiently many particles, specifically with N much greater than O(R0U/ν), the initial evolution is strongly influenced by wake-mediated interactions; here, N is the total number of particles, and R0 is the initial cluster radius. As a result, the cluster first evolves into a nearly planar configuration with an asymptotically small aspect ratio of O(R0U/N ν), the plane of the cluster being perpendicular to the direction of gravity; subsequent expansion occurs with an unchanged aspect ratio. For relatively sparse clusters with N smaller than O(R0U/ν), the probability of wake interactions remains negligible, and the cluster expands while retaining its spherical shape. The long-time expansion in the former case, and that for all times in the latter case, is driven by disturbance velocity fields produced by the particles outside their wakes. The resulting interactions between particles are therefore mutually repulsive with forces that obey an inverse-square law. The analysis presented describes cluster evolution in this regime. A continuum representation is adopted with the clusters being characterized by a number density field (n(r, t)), and a corresponding induced velocity field (u (r, t)) arising on account of interactions. For both planar axisymmetric clusters and spherical clusters with radial symmetry, the evolution equation admits a similarity solution; either cluster expands self-similarly for long times. The number density profiles at different times are functions of a similarity variable η = (r/t1/3), r being the radial distance away from the cluster centre, and t the time. The radius of the expanding cluster is found to be of the form Rcl (t) = A (ν a)1/3N1/3t1/3, where the constant of proportionality, A, is determined from an analytical solution of the evolution equation; one finds A = 1.743 and 1.651 for planar and spherical clusters, respectively. The number density profile in a planar axisymmetric cluster is also obtained numerically as a solution of the initial value problem for a canonical (Gaussian) initial condition. The numerical results compare well with theoretical predictions, and demonstrate the asymptotic stability of the similarity solution in two dimensions for long times, at least for axisymmetric initial conditions.