  Phoretic particles exploit local self-generated physico-chemical gradients to
achieve self-propulsion at the micron scale. The collective dynamics of a large
number of such particles is currently the focus of intense research efforts,
both from a physical perspective to understand the precise mechanisms of the
interactions and their respective roles, as well as from an experimental point
of view to explain the observations of complex dynamics as well as formation of
coherent large-scale structures. However, an exact modelling of such
multi-particle problems is difficult and most efforts so far rely on the
superposition of far-field approximations for each particle's signature, which
are only valid asymptotically in the dilute suspension limit. A systematic and
unified analytical framework based on the classical Method of Reflections (MoR)
is developed here for both Laplace and Stokes' problems to obtain the
higher-order interactions and the resulting velocities of multiple phoretic
particles, up to any order of accuracy in the radius-to-distance ratio
$\varepsilon$ of the particles. Beyond simple pairwise chemical or hydrodynamic
interactions, this model allows us to account for the generic
chemo-hydrodynamic couplings as well as $N$-particle interactions ($N\geq 3$).
The $\varepsilon^5$-accurate interaction velocities are then explicitly
obtained and the resulting implementation of this MoR model is discussed and
validated quantitatively against exact solutions of a few canonical problems.
