Integral with messy integrand Let $N\in\mathbb{N}$, $\theta_{i}>0$, and $a_{ij}\in\mathbb{R}$, $\forall i = 1,\ldots, N$. Is there a somewhat explicit expression for
\begin{align}
\int_{0}^{\infty} t \left[\prod_{i=1}^{N} \frac{1}{\sqrt{1+2\theta_{i}t}}\right] 
    \left[\sum_{i,j=1}^{N}\frac{a_{ij}}{(1+2\theta_{i}t)(1+2\theta_{j}t)}\right]dt?
\end{align}
 A: When $N>2$, assuming that all $\theta_i$'s are different, we can derive an expression for the integral which involves a sum of multivariate hypergeometric functions.
We decompose the terms of the summation, considering first the non-diagonal ones:
\begin{align}
 \sum_{\stackrel{i,j=1}{i\ne j}}^{N}\frac{a_{ij}}{(1+2\theta_{i}t)(1+2\theta_{j}t)}&=\sum_{\stackrel{i,j=1}{i\ne j}}^{N}\frac{a_{ij}}{\theta_i-\theta_j}\left[\frac{\theta_i}{1+2\theta_{i}t}-\frac{\theta_j}{1+2\theta_{j}t}\right]\\
 &=\sum_{\stackrel{i,j=1}{i\ne j}}^{N}\left( a_{ij}+a_{ji} \right)\frac{\theta_i}{\theta_i-\theta_j}\frac{1}{1+2\theta_{i}t}\\
 &=\sum_{i=1}^N \frac{\alpha_i}{1+2\theta_{i}t}
\end{align}
where we define
\begin{equation}
 \alpha_i=\sum_{\stackrel{j=1}{j\ne i}}^{N}\left( a_{ij}+a_{ji} \right)\frac{\theta_i}{\theta_i-\theta_j}
\end{equation} 
(all $\theta_i$ are supposed to be different). Taking into account the diagonal contribution, the summation becomes
\begin{align}
\sum_{i,j=1}^{N}\frac{a_{ij}}{(1+2\theta_{i}t)(1+2\theta_{j}t)}&= \sum_{i=1}^N \left[\frac{a_{ii}}{\left( 1+2\theta_{i}t \right)^2}+\frac{\alpha_i}{1+2\theta_{i}t}\right]\\
\end{align} 
Then
\begin{align}
  \prod_{k=1}^{N} &\frac{1}{\sqrt{1+2\theta_{k}t}}\sum_{i,j=1}^{N}\frac{a_{ij}}{(1+2\theta_{i}t)(1+2\theta_{j}t)}\\
  &=
  \sum_{i=1}^N \left[\frac{a_{ii}}{\left( 1+2\theta_{i}t \right)^{5/2}}+\frac{\alpha_i}{\left( 1+2\theta_{i}t \right)^{3/2}}\right]
  \prod_{\stackrel{k=1}{k\ne i}}^{N} \left( 1+2\theta_{k}t \right)^{-1/2}\\
\end{align} 
A multivariate hypergeometric function is defined as
\begin{align}
 R_{-a}\left( \mathbf{b};\mathbf{z} \right)&=R_{-a}\left(b_1,b_2,\ldots,b_n;z_1,z_2,\ldots,z_n\right)\\
 &=\frac{1}{\mathrm{B}\left(a,a^{\prime}\right)}\int_{0}^{\infty}t^{a%
-1}\prod^{n}_{j=1}(1+tz_{j})^{-b_{j}}\mathrm{d}t
\end{align} 
where
 $b_1+b_2+\cdots+b_n>a>0, b_j\in\mathbb{R},z_j\in\mathbb{C}\backslash (-\infty,0]$, and
\begin{equation}
 a^{\prime}=-a+\sum_{j=1}^{n}b_{j}
\end{equation}
The integral can be written as,
\begin{align}
 I&=\int_{0}^{\infty} t \left[\prod_{k=1}^{N} \frac{1}{\sqrt{1+2\theta_{k}t}}\right] 
    \left[\sum_{i,j=1}^{N}\frac{a_{ij}}{(1+2\theta_{i}t)(1+2\theta_{j}t)}\right]dt\\
    &=\frac{1}{4}\sum_{i=1}^N \int_{0}^{\infty} t\left[\frac{a_{ii}}{\left( 1+\theta_{i}t \right)^{5/2}}+\frac{\alpha_i}{\left( 1+\theta_{i}t \right)^{3/2}}\right]
  \prod_{\stackrel{k=1}{k\ne i}}^{N} \left( 1+\theta_{k}t \right)^{-1/2}\,dt\\
  &=\frac{1}{N}\sum_{i=1}^N\left[\frac{a_{ii}}{N+2}R_{-2}\left( \mathbf{\frac{1}{2}}+2\mathbf{e_i};\boldsymbol{\theta} \right)+\frac{\alpha_{i}}{N-2}R_{-2}\left( \mathbf{\frac{1}{2}}+\mathbf{e_i};\boldsymbol{\theta} \right)\right]
\end{align} 
Here $\mathbf{e_i}$ is an n-tuple with 1 in the i-th place and 0's elsewhere.
The multivariate hypergeometric function can be expressed in terms of a Lauricella $F_D$ function, i.e. as a hypergeometric series of the variables $1-\theta_i$ (see for example this paper by B.C. Carlson).
When several $\theta_i$ are identical, the above method must be modified by taking into account a lower number of variables and different exponent.  Finally, the cases $N=1,2$ can be calculated directly, they lead to elliptic integral expressions. 
