The cumulative distribution function is
$$P(Y_1 \leq y_1, Y_2 \leq y_2) = H_{\alpha}(y_1,y_2) =
\log_{\alpha} [1 + (\alpha^{y_1}-1)(\alpha^{y_2}-1)/
(\alpha-1)] $$
for $alpha != 1$.
Note the logarithm here is to base $alpha$.
The support of the function is the unit square.
When $01$ then
$h_{1/alpha}(1-y_1,y_2)$.
If $alpha=1$ then $H(y1,y2)=y1*y2$,
i.e., uniform on the unit square.
As $alpha$ approaches 0 then
$H(y1,y2)=min(y1,y2)$.
As $alpha$ approaches infinity then
$H(y1,y2)=max(0,y1+y2-1)$.
The default is to use Fisher scoring implemented using
rbifrankcop.
For intercept-only models an alternative is to set nsimEIM=NULL
so that a variant of Newton-Raphson is used.