The binding of a system of N polarons subject to a constant magnetic field of strength B is investigated within the Pekar–Tomasevich approximation. In this approximation, the energy of N polarons is described in terms of a non-quadratic functional with a quartic term that accounts for the electron–electron self-interaction mediated by phonons. The size of a coupling constant, denoted by α, in front of the quartic term is determined by the electronic properties of the crystal under consideration, but in any case it is constrained by 0 < α < 1. For all values of N and B, we find an interval α N,B < α < 1 where the N polarons bind in a single cluster described by a minimizer of the Pekar–Tomasevich functional. This minimizer is exponentially localized in the N-particle configuration space $${\mathbb{R}^{3N}}$$ .