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
rbifrankcop.
For intercept-only models an alternative is to set nsimEIM=NULL
so that a variant of Newton-Raphson is used.