  The recent gravitational wave transient GW190521 has been interpreted by the
LIGO-Virgo collaboration (LVC) as sourced by a binary black hole (BH) merger.
According to the LVC parameter estimation, at least one of these progenitors
falls into the so-called pair-instability supernova mass gap. This raises the
important question of how and when these progenitors formed. In this paper we
use an accretion model with super-Eddington mass accretion rate obtained from
General Relativity hydrodynamics simulations to analyse the scenario wherein
the GW190521 original progenitors (OPs) formed at lower masses (and spins) and
grew to their estimated LVC parameters by relativistic accretion. We consider
that the environment wherein the binary is immersed has density gradients as
well as a dependence on the Mach number of the gas. Taking the LVC parameter
estimation at $z=0.82$ as the endpoint of the accretion evolution, we estimate
the initial masses and spins of the OPs at three different red-shifts $z=100, \
50$, and $20$. We found three distinct possible types of OPs: $(i)$ $10^{-4}
M_{\odot} - 3 M_{\odot}$ almost non-rotating (with Kerr spin parameter
$a_{\star}< 10^{-2}$) primordial BHs; $(ii)$ $3 M_{\odot} - 40M_{\odot}$ slowly
rotating ($ 10^{-2} < a_{\star} < 0.5$) stellar mass BHs; $(iii)$ $40M_{\odot}
- 70M_{\odot}$ BHs with a moderate spin parameter $a_{\star}\sim 0.5$, which
could originate from the collapse of high mass Pop III stars. The mass spread
is due to varying the density gradient and the relativistic Mach number of the
cosmic plasma; the variation of the masses due to the origin at different
red-shifts, on the other hand, is negligible, $\sim 2\%$ ...
