Archimedean Clayton copula entropy Question
I would like to derive the entropy of Archimedean parametric copulas  (Clayton, Frank, or Gumbel), focusing here on the Clayton copula.
Link to similar question on the t-copula.
The bivariate copula function, $C$, for the Clayton copula, with transformed marginals $u$ and $v$, and dependence parameter $\theta\in \mathbb{R}_{\geq 0}$, is
$$ C(u, v) = \bigg[ u^{-\theta} + v^{-\theta} -1 \bigg]^{-1/\theta} $$
Its copula density is the second mixed partial derivative of $C(u,v)$:
\begin{align}
c\left(u,v ; \theta\right) = (1+\theta)(u \cdot v)^{-1-\theta}(u^{-\theta}+v^{-\theta}-1)^{-\frac{1}{\theta}-2}
\end{align}
The differential entropy of a univariate variable's pdf, $f(x)$, is typically
$$H(X)=-\int_{-\infty} ^{\infty} f(x) \ln f(x) dx$$
while any copula entropy can be estimated as
$$H(c(u,v))=-\int_{[0,1]^2} c(u,v) \ln c(u,v) \hspace{1mm} du \hspace{1mm} dv $$
How can we derive a closed-form analytical solution for the entropy of the Clayton copula density $c(u,v)$?
First attempt
\begin{align}
H(c(u,v)) =&  -\iint_{[0,1]^2} c\left(u, v ; \theta\right) \ln c\left(u, v ; \theta\right) \mathrm{d}u \, \mathrm{d}v \\
= &  -\iint_{[0,1]^2} (1+\theta)\left(u v\right)^{-1-\theta}\left(u^{-\theta}+v^{-\theta}-1\right)^{-\frac{1}{\theta}-2}\\
\times & \ln\left[(1+\theta)\left(u v\right)^{-1-\theta}\left(u^{-\theta}+v^{-\theta}-1\right)^{-\frac{1}{\theta}-2}\right] \mathrm{d}u \, \mathrm{d}v \\
= &  -(1+\theta) \iint_{[0,1]^2} \left(u v\right)^{-1-\theta}\left(u^{-\theta}+v^{-\theta}-1\right)^{-\frac{1}{\theta}-2}\\
\times & \left[\ln(1+\theta) - (1+\theta)\ln\left(u v\right) -\left(\frac{1}{\theta}+2\right)\ln\left(u^{-\theta}+v^{-\theta}-1\right)\right] \mathrm{d}u \, \mathrm{d}v \\
= &  -(1+\theta) \iint_{[0,1]^2} \bigg[\ln(1+\theta) c - (1+\theta)\ln\left(u v\right) c -\left(\frac{1}{\theta}+2\right)\ln\left(u^{-\theta}+v^{-\theta}-1\right) c \bigg] \ \mathrm{d}u \, \mathrm{d}v \\
= &  -(1+\theta) \ln(1+\theta) \iint_{[0,1]^2}  c \ \mathrm{d}u \, \mathrm{d}v + (1+\theta)^2\iint_{[0,1]^2} c \ln\left(u v\right)  \ \mathrm{d}u \, \mathrm{d}v \\
+ & (1+\theta) \left(\frac{1}{\theta}+2\right) \iint_{[0,1]^2} c \ln\left(u^{-\theta}+v^{-\theta}-1\right) \ \mathrm{d}u \, \mathrm{d}v \\
\stackrel{\dagger}{=} &  -(1+\theta) \ln(1+\theta) \int_0^1  c  \bigg|_{v=0}^{v=1} \ \mathrm{d}u + (1+\theta)^2 \int_0^1 c \ln\left(u v\right) \bigg|_{v=0}^{v=1} \ \mathrm{d}u\\
+ & (1+\theta) \left(\frac{1}{\theta}+2\right) \int_0^1 c \ln\left(u^{-\theta}+v^{-\theta}-1\right) \bigg|_{v=0}^{v=1} \ \mathrm{d}u  \\
= & ?
\end{align}
where $c = \left(u v\right)^{-1-\theta}\left(u^{-\theta}+v^{-\theta}-1\right)^{-\frac{1}{\theta}-2} $, and $\dagger$ is due to treating the double integrals as an iterated integral,
$$\iint_{[0,1]^2}  c(u,v) \ \mathrm{d}u \, \mathrm{d}v  = \int_0^1  c(u,v)   \bigg|_{v=0}^{v=1} \mathrm{d}u $$
Why do I think an analytical solution of copula entropy can be found? Because there is one for the entropy of the Normal distribution's pdf, derived here.
 A: The place for a comment was too short for the following, so this became an answer.
Well, removing a factor (namely the factor that makes live too short to compute...) simplifies things in a high measure, but we still have a mess... For the above i see only some first steps, but then things still get complicated. The $\theta$ is hard to type, there will be a $t$ instead.
We have
$$
\begin{aligned}
H
&=-
\iint_{[0,1]^2} 
\frac
{(1+t)(uv)^{-t}}
{(u^{-t} + v^{-t} -1)^{2+1/t}}\;\ln 
\frac
{(1+t)(uv)^{-t-1}}
{(u^{-t} + v^{-t} -1)^{2+1/t}}
\;
\frac {du}u \; \frac {dv}v
\\
&= 
\iint_{I^2} 
\frac
{(1+t)UV}
{(U + V -1)^{2+1/t}}
\;
\ln 
\frac
{(1+t)(UV)^{-(1+t)/t}}
{(U + V-1)^{2+1/t}}
\;
\frac 1{t^2}
\;
\frac {dU}U \; \frac {dV}V
\\
&= 
\frac {1+t}{t^2}
\iint_{I^2} 
\frac
1{(U + V -1)^{2+1/t}}
\;\ln 
\frac
{(1+t)(UV)^{-(1+t)/t}}
{(U + V - 1)^{2+1/t}}
\;
dU\; dV\ .
\end{aligned}
$$
We substituted $U=u^{-t}$, $V=v^{-t}$, so $\frac {dU}U$ is
$(-t)\frac {du}u$, and $\frac {dV}V$ is $(-t)\frac {dv}v$, so that we obtain a better looking expression.
The integral is now over $I^2$, where $I$ is $[1,\infty]$, because of the sign of $-t$ in $U=u^{-t}$.
Now under the logarithm we split four factors. And have to compute correspondingly four integrals.

*

*The integral in $\ln(t+1)$ is the simplest. It is in fact an integral in $W=(U+V)\ge 2$. We have working with $a=2+1/t$
$$\iint_{I^2} 
\frac
1{(U + V -1)^a}
\; dU\; dV
=
\int_2^\infty\frac{W-2}{(W-1)^a}\; dW
\\
=\frac 1{a^2-3a+2}
=\frac 1{a-2}-\frac 1{a-1}
=
t-\frac t{1+t}
\ .
$$

*The integrals in $\ln U$ and $\ln V$ are equal, it is enough to compute only one of them. Again pass from $(U,V)$ to $(U,W)$, where $W=U+V$ formally. So we have to integrate something like
$$
\iint_{\substack{1\le U<\infty\\2\le W<\infty\\1+U\le W}}
\frac
1{(W -1)^a}
\ln U
\; dU\; dW
$$
We can first integrate in $U$ from $1$ to $W-1$, but some $\log(W-2)$ will come into play and the work begins.

*The final integral can also be arranged in a better form, passing from $(U,V)$ to $(U,W)$ as above, but then we have to compute something like
$$
\int_{2}^\infty
\frac{W-2}{(W-1)^a}\ln(W-1)\;dW\ .
$$

Now try to apply integration by parts to get rid of the logarithmic term.
For special values of $t$ (and $a$) this may be computed, but i will stop here.
