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 \ne 1$.
Note the logarithm here is to base $\alpha$.
The support of the function is the unit square. When $0 < \alpha < 1$ the probability density function
$h_{\alpha}(y_1,y_2)$
is symmetric with respect to the lines $y_2=y_1$
and $y_2=1-y_1$.
When $\alpha > 1$ then
$h_{\alpha}(y_1,y_2) = h_{1/\alpha}(1-y_1,y_2)$.
If $\alpha=1$ then $H(y_1,y_2) = y_1 y_2$,
i.e., uniform on the unit square.
As $\alpha$ approaches 0 then
$H(y_1,y_2) = \min(y_1,y_2)$.
As $\alpha$ approaches infinity then
$H(y_1,y_2) = \max(0, y_1+y_2-1)$.
The default is to use Fisher scoring implemented using
rfrank.
For intercept-only models an alternative is to set nsimEIM=NULL
so that a variant of Newton-Raphson is used.