  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 ($m_1$) distribution are consistent
with the gravitational wave (GW) observations. Notably, Pop III and EMP ($<
10^{-2}$ $Z_\odot$) binary stars yield most of the pair instability (PI) mass
gap events with $m_1 = 65$--$130$ $M_\odot$. Pop III binary stars contribute
more to the PI mass gap events with increasing redshift, and all the PI mass
gap events have the Pop III origin at redshifts $\gtrsim 8$. Our result can be
assessed by future GW observations in the following two points. First, there
are no binary BHs with $m_1=100$--$130$ $M_\odot$ in our result, and thus the
$m_1$ distribution should suddenly drop in the range of $m_1=100$--$130$
$M_\odot$. Second, the PI mass gap event rate should increase toward higher
redshift up to $\sim 11$, since those events mainly originate from the Pop 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 Pop 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.
