The detection of double black hole (BH+BH) mergers provides a unique possibility to understand their physical properties and origin. To date, the LIGO–Virgo–KAGRA network of high-frequency gravitational-wave observatories has announced the detection of more than 85 BH+BH merger events. An important diagnostic feature that can be extracted from the data is the distribution of effective inspiral spins of the BHs. This distribution is in clear tension with theoretical expectations from both an isolated binary star origin, which traditionally predicts close-to-aligned BH component spins, and formation via dynamical interactions in dense stellar environments that predicts a symmetric distribution of effective inspiral spins. Here it is demonstrated that isolated binary evolution can convincingly explain the observed data if BHs have their spin axis tossed during their formation process in the core collapse of a massive star, similarly to the process evidently acting in newborn neutron stars. BH formation without spin-axis tossing, however, has difficulties reproducing the observed data—even if alignment of spins prior to the second core collapse is disregarded. Based on simulations with only a minimum of assumptions, constraints from empirical data can be made on the spin magnitudes of the first- and second-born BHs, thereby serving to better understand massive binary star evolution prior to the formation of BHs.