  Our understanding of the formation and evolution of binary black holes (BBHs)
is significantly impacted by the recent discoveries made by the LIGO/Virgo
collaboration. Of utmost importance is the detection of the most massive BBH
system, GW190521. Here we investigate what it takes for field massive stellar
binaries to account for the formation of such massive BBHs. Whether the high
mass end of the BH mass function is populated by remnants of massive stars that
either formed at extremely low metallicities and avoid the pair-instability
mass gap or increase their birth mass beyond the pair-instability mass gap
through the accretion of gas from the surrounding medium. We show that assuming
that massive stars at very low metallicities can form massive BHs by avoiding
pair-instability supernova, coupled with a correspondingly high formation
efficiency for BBHs, can explain the observed BH mass function. To this end,
one requires a relation between the initial and final mass of the progenitor
stars at low metallicities that is shallower than what is expected from wind
mass loss alone. On the other hand, assuming pair-instability operates at all
metallicities, one can account for the observed BH mass function if at least
about 10% of the BHs born at very low metallicities double their mass before
they merge because of accretion of ambient gas. Such BBHs will have to spend
about a Gyr within a parsec length-scale of their parent atomic cooling halos
or a shorter timescale if they reside in the inner sub-parsecs of their host
dark matter halos. Future stellar evolution calculations of massive stars at
very low metallicity and hydrodynamical simulations of gas accretion onto BBHs
born in atomic cooling halos can shed light on this debate.
