ABSTRACT The existence of massive stellar black hole binaries (MBHBs), with primary black hole masses $\ge 31 \, \mathrm{ M}_⊙$, was proven by the detection of the gravitational wave (GW) event GW150914 during the first LIGO/Virgo observing run (O1), and successively confirmed by seven additional GW signals discovered in the O1 and O2 data. By adopting the galaxy formation model gamesh coupled with binary population synthesis (BPS) calculations, here we investigate the origin of these MBHBs by selecting simulated binaries compatible in mass and coalescence redshifts. We find that their cosmic birth rates peak in the redshift range 6.5 ≤ z ≤ 10, regardless of the adopted BPS. These MBHBs are then old systems forming in low-metallicity ($Z ∼ [0.01\!-\!0.1] \, Z_{⊙ }$), low-stellar-mass galaxies, before the end of cosmic reionization, i.e. significantly beyond the peak of cosmic star formation. GW signals generated by coalescing MBHBs open up new possibilities to probe the nature of stellar populations in remote galaxies, at present too faint to be detected by available electromagnetic facilities.