Definite integral closed-form expression

Integrating by parts we have:

F[t_] := Exp[-a*t];
G[t_] := Log[t]*Log[1 + t];
HoldForm[Integrate[F[t]*G'[t], t] = F[t]*G[t] - Integrate[F'[t]*G[t], t]]
integral = Integrate[F[t]*G'[t], t] == F[t]*G[t] - Integrate[F'[t]*G[t], t] 
First@Expand@Solve[integral, \[Integral]E^(-a t) Log[t] Log[1 + t] \[DifferentialD]t]

$$\int F(t) G'(t) \, dt=F(t) G(t)-\int F'(t) G(t) \, dt$$ $$\int e^{-a t} \left(\frac{\log (t)}{t+1}+\frac{\log (t+1)}{t}\right) \, dt=e^{-a t} \log (t) \log (t+1)+a \int e^{-a t} \log (t) \log (t+1) \, dt$$ $$\int e^{-a t} \log (t) \log (t+1) \, dt= \frac{\int e^{-a t} \left(\frac{\log (t)}{t+1}+\frac{\log (t+1)}{t}\right) \, dt}{a}-\frac{e^{-a t} \log (t) \log (t+1)}{a}$$

define integral is then:

Limit[-Exp[-a*t]*Log[t]*Log[1 + t]/a, t -> Infinity, Assumptions -> a > 0] - 
Limit[-Exp[-a*t]*Log[t]*Log[1 + t]/a, t -> 0, Assumptions -> a > 0] +
Integrate[Log[1 + t]/(Exp[a*t]*t), {t, 0, Infinity},Assumptions -> a > 0]/a + 
Integrate[Log[t]/(Exp[a*t]*(1 + t)), {t, 0, Infinity}, Assumptions -> a > 0]/a

$$\frac{G_{2,3}^{3,1}\left(a\left| \begin{array}{c} 0,1 \\ 0,0,0 \\ \end{array} \right.\right)-e^a \left(G_{2,3}^{3,0}\left(a\left| \begin{array}{c} 1,1 \\ 0,0,0 \\ \end{array} \right.\right)+(\log (a)+\gamma ) \Gamma (0,a)\right)}{a}$$

Numerical check for: a=2

a = 2;
NIntegrate[Exp[-a*t]*(Log[t]*Log[t + 1]), {t, 0, Infinity}]
(*-0.0730318*)

1/a (-E^a (Gamma[0, a] (EulerGamma + Log[a]) + MeijerG[{{}, {1, 1}}, {{0, 0, 0}, {}}, a]) + MeijerG[{{0}, {1}}, {{0, 0, 0}, {}}, a])//N
(*-0.0730318*)

EDIT

I have now confirmed my closed form expression numerically. This was possible by helping Mathematica to calculate the numerical value of mixed partial derivatives of HypergeometricU[a,b,z] with respect to a and b.

Original post

We derive a closed form by another procedure. The procedure seems to be valid but Mathematica has difficulties with the numeric values of the derivatives of HypergeometricU[a,b,z] with respect to a and b.

Nevertheless it might be useful to have a look at these short developments the results of which resemble those of Mariusz.

$Version

(* Out[7]= "10.1.0  for Microsoft Windows (64-bit) (March 24, 2015)" *)

Consider the related integral

f = 
 Integrate[Exp[-a t ] t^p (1 + t)^q, {t, 0, \[Infinity]}, 
  Assumptions -> {a > 0, p > -1}]

(* Out[21]= Gamma[1 + p] HypergeometricU[1 + p, 2 + p + q, a] *)

Now we retrieve the original integral by twofold differentiation

D[f,p]/.p->0

$$U^{(0,1,0)}(1,q+2,a)+U^{(1,0,0)}(1,q+2,a)+\gamma e^a \left(-a^{-q-1}\right) \Gamma (q+1,a)$$

D[%, q] /. q -> 0

$$-\frac{\gamma e^a \left(G_{1,2}^{2,0}\left(a\left| \begin{array}{c} 1 \\ 0,0 \\ \end{array} \right.\right)+e^{-a} \log (a)\right)}{a}+U^{(0,2,0)}(1,2,a)+U^{(1,1,0)}(1,2,a)+\frac{\gamma \log (a)}{a}$$

FunctionExpand[%]

$$U^{(0,2,0)}(1,2,a)+U^{(1,1,0)}(1,2,a)-\frac{\gamma e^a \left(-\text{Ei}(-a)+\frac{1}{2} \left(\log (-a)-\log \left(-\frac{1}{a}\right)\right)+e^{-a} \log (a)-\log (a)\right)}{a}+\frac{\gamma \log (a)}{a}$$

Following the comment of Alex we finally impose a>0 and find

ff1= Simplify[ff, a > 0]

$$U^{(0,2,0)}(1,2,a)+U^{(1,1,0)}(1,2,a)+\frac{\gamma e^a \text{Ei}(-a)}{a}$$

This is the closed form expression of the integral in terms of special functions as requested in the OP.

Numerical check (part of EDIT)

We now check the numerical agreement. As mentioned before Mathematica refuses to evaluate the expression.

In order to see in more detail which part makes the problem we split the expression into a list

ff2 = List @@ ff1

$$\left\{\frac{\gamma e^a \text{Ei}(-a)}{a},U^{(0,2,0)}(1,2,a),U^{(1,1,0)}(1,2,a)\right\}$$

 % /. a -> 1.

$$\left\{-0.344221,0.531931,U^{(1,1,0)}(1,2,1.)\right\}$$

As Mathematica does not evaluate the mixed derivative numerically we give it a little help defining the partial derivative explicitly as a quotient of diffences in a sufficiently small interval

ff13[a_] := ff13[a_] := 
 With[{d = 10^(-3)}, (1/
      d)*(Derivative[1, 0, 0][HypergeometricU][p, q + d/2, a] - 
          Derivative[1, 0, 0][HypergeometricU][p, q - d/2, a]) /. {p -> 
     1, q -> 2}]

Where smaller d does not change the result.

The explicit expression of the integral which can also be evaluated numerically is given by

ff1N[a_] = ff1[[1]] + ff1[[2]] + ff13[a];

Taking the same value as Mariusz did we find

ff1N[2.]

(* Out[191]= -0.0730318 *)

For comparison the numerical integral is

fi[a_] := NIntegrate[Exp[-a t] Log[t] Log[1 + t], {t, 0, \[Infinity]}]

There is excellent agreement over a greater range of "a" as can be checked with this plot (to be done by the reader as it just shows a line y = 1)

Plot[{fi[a] ff1N[a]}, {a, 0, 2}, PlotRange -> {0.99, 1.01}] (* plot not shown *)]