A boundary-collocation technique, used earlier by Chen and Keh [(1996) Aerosol Sci. Technol. 24, 21-35] in the study of the thermophoresis of N coaxial spheres along the line of their centers, is extended to describe the motion of an assemblage of N aerosol spheres arranged arbitrarily in three-dimensional space. The spheres are allowed to differ in radius, in thermal conductivity and in surface properties; they may move independently, or they may be linked by infinitesimally thin rods to form a rigid aggregate. The Knudsen numbers are assumed to be small so that the fluid flow is described by a continuum model with a thermal creep and a hydrodynamic slip at the particle surfaces. Results are presented in terms of pair-interaction coefficients for the thermophoretic velocities of the particles. For two-sphere systems, the translational and angular velocities of the particles at all orientations and separation distances agree very well with the exact solutions obtained by using spherical bipolar coordinates or the asymptotic solutions obtained by using a method of reflections. The particle-interaction parameters of linear chains of three spheres show that the existence of the third sphere can significantly affect the mobilities of the other two spheres. For the cases of a rigid dumbbell composed of two spheres, the numerical solutions for the particle velocities compare quite favorably with the formulas derived analytically. Finally, our numerical results for the interaction between two spheres are used to find the effect of volume fractions of particles of each type on the mean thermophoretic velocities in a polydisperse aerosol.