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}

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.