  Using the 3D smoothed particle hydrodynamics code PHANTOM, we investigate the
evolution of the orbital properties of massive black hole binaries embedded in
massive discs where gravitational instabilities (GIs) triggered by the disc
self-gravity are the only source of angular momentum transport. In particular,
we investigate the evolution of binaries with different initial eccentricities
$e_0=0.05,\,0.5,\,0.8$ and mass ratios $q=0.1,\,0.3,\,0.9$. Our simulations
suggest that there might not be a unique value of critical eccentricity. We
find initially more eccentric binaries to tend to higher asymptotic
eccentricity values than more circular ones. This implies that there is a range
of critical eccentricity values, that depends on the initial condition of the
system. In particular, we find the width of this range to be narrower for more
unequal binaries. We furthermore measure preferential accretion onto our
binaries, finding more accretion onto the primary only for mass ratio $q=0.3$
and eccentricity $e=0.8$. We discuss how this might have implications for the
amplitude of the gravitational wave background detected by Pulsar Timing Arrays
(PTA) experiments. We finally measure the corresponding value of the viscosity
parameter $\alpha$ in our simulations and discuss how this depends on the
binary properties.
