<title>Abstract</title>
We investigate the formation of merging binary black holes (BHs) through isolated binary evolution, performing binary population synthesis calculations covering an unprecedentedly wide metallicity range of Population (Pop) I, II, III, and extremely metal-poor (EMP) binary stars. We find that the predicted merger rate density and primary BH mass (<italic>m</italic>
₁) distribution are consistent with the gravitational wave (GW) observations. Notably, Population III and EMP (&lt;10⁻²
<italic>Z</italic>
_⊙) binary stars yield most of the pair instability (PI) mass gap events with <italic>m</italic>
₁ = 65–130 <italic>M</italic>
_⊙. Population III binary stars contribute more to the PI mass gap events with increasing redshift, and all the PI mass gap events have the Population III origin at redshifts ≳8. Our result can be assessed by future GW observations in the following two points. First, there are no binary BHs with <italic>m</italic>
₁ = 100–130 <italic>M</italic>
_⊙ in our result, and thus the <italic>m</italic>
₁ distribution should suddenly drop in the range of <italic>m</italic>
₁ = 100–130 <italic>M</italic>
_⊙. Second, the PI mass gap event rate should increase toward higher redshift up to ∼11, since those events mainly originate from the Population III binary stars. We find that the following three assumptions are needed to reproduce the current GW observations: a top-heavy stellar initial mass function and the presence of close binary stars for Population III and EMP binary stars, and inefficient convective overshoot in the main-sequence phase of stellar evolution. Without any of the above, the number of PI mass gap events becomes too low to reproduce current GW observations.