Difficult Laplace Transform Type of Integral

I have only a partial answer: it pretty much just consists of repeated integration by parts, so maybe this isn't what you're looking for.

Write $A(n_1, \dots, n_r) = \int_0^\infty e^{-as} \prod_i T_{n_i}(s) \,ds$. Then integration by parts gives \begin{align*} A(n_1, \dots, n_r) &= -\frac{1}{a}\bigg[e^{-as}\prod_i T_{n_i}(s)\bigg]\bigg|_0^\infty + \frac{1}{a}\sum_{j : n_j > 0} \int_0^\infty e^{-as} T_{n_j - 1}(s)\prod_{i \neq j} T_{n_i}(s) \,ds \\ &= \frac{1}{a} + \frac{1}{a} \sum_{j : n_j > 0} A(n_1, \dots, n_j - 1, \dots, n_r) \end{align*} when all $n_i \geq 0$.

To make things a little easier, let's take $T_n(s) = 0$ when $n < 0$ (this is consistent with the identity $T_n'(s) = T_{n-1}(s)$), so we'll have $A(n_1, \dots, n_r) = 0$ when some $n_i < 0$, and so we can write $$A(n_1, \dots, n_r) = \frac{1}{a}\bigg(\mathbb{1}_{(n_i)} + \sum_i A(n_1, \dots, n_i - 1, \dots, n_r)\bigg)$$ for any integers $n_1, \dots, n_r$, where $\mathbb{1}_{(n_i)} = 1$ if all $n_i \geq 0$ and is $0$ otherwise.

Letting $f(x_1, \dots, x_r) = \sum A(n_1, \dots, n_r) x_1^{n_1} \cdots x_r^{n_r}$ be the generating function of $A$, this identity gives $$af = \frac{1}{(1-x_1) \cdots (1-x_r)} + (x_1 + \cdots + x_r)f$$ hence \begin{align*} f &= \frac{1}{(1 - x_1) \cdots (1 - x_r)(a - x_1 - \cdots - x_r)} \\ &= \frac{1}{a} \sum_{l_1, \dots, l_r, m_1, \dots, m_r} a^{-m_1 - \cdots - m_r}\binom{m_1 + \cdots + m_r}{m_1, \dots, m_r} x_1^{l_1 + m_1} \cdots x_r^{l_r + m_r} \\ &= \frac{1}{a} \sum_{n_1, \dots, n_r} \sum_{0 \leq m_i \leq n_i} a^{-m_1 - \cdots - m_r}\binom{m_1 + \cdots + m_r}{m_1, \dots, m_r} x_1^{n_1} \cdots x_r^{n_r} \\ \end{align*} and finally, $$A(n_1, \dots, n_r) = \frac{1}{a} \sum_{0 \leq m_i \leq n_i} a^{-m_1 - \cdots - m_r}\binom{m_1 + \cdots + m_r}{m_1, \dots, m_r} $$ (this could alternatively be found by repeatedly expanding the recurrence).

The case relevant to the question is the one where $n_1 = \cdots = n_r = N$, $r = k-j$, $a = 1/x - j$, but I'm not sure how to simplify this sum. I've been thinking about it in terms of lattice paths in $r$ dimensions: in particular the coefficient of $a^{-t}$ is the number of lattice paths of length $t$ where we have at most $n_i$ steps in the $i$-th direction for each $i$, but that doesn't seem to help simplify any further.

As noted by DinosaurEgg, we're only interested in $$A(N, \dots, N) = \frac{1}{a} \sum_{0 \leq m_i \leq N} a^{-m_1 - \cdots - m_r}\binom{m_1 + \cdots + m_r}{m_1, \dots, m_r}$$ but this sum still seems difficult to simplify to me, because of the $m_i \leq N$ constraints.

So long as $j < 1/x$ it is possible to write the function as a finite sum of terms (i.e., in closed form) by writing the power of the truncated exponential as a polynomial, and then turning the integral into a polynomial sum of gamma functions. If we apply the multinomial theorem we can write the power of the truncated exponential as a polynomial of degree $N(k-j)$ as follow:

$$\begin{equation} \begin{aligned} \left[ \sum_{n=0}^N \frac{s^n}{n!} \right]^{k-j} &= \sum_{r_0 + \cdots + r_N = k-j} {k-j \choose \mathbf{r}} \prod_{n=0}^N \Big( \frac{s^n}{n!} \Big)^{r_n} \\[6pt] &= \sum_{r_0 + \cdots + r_N = k-j} \frac{(k-j)!}{\prod_{n=0}^N (r_n!) (n!)^{r_n}} \prod_{n=0}^N s^{n \cdot r_n} \\[6pt] &= \sum_{r_0 + \cdots + r_N = k-j} \frac{(k-j)!}{\prod_{n=0}^N (r_n!) (n!)^{r_n}} s^{\sum_{n=0}^N n \cdot r_n} \\[6pt] &= \sum_{i=0}^{N(k-j)} a_i s^i, \\[6pt] \end{aligned} \end{equation}$$

where the coefficients $a_0,a_1,...,a_{N(k-j)}$ depend on $N$ and $k-j$, and are obtained from this multinomial summation. We then have:

$$\begin{equation} \begin{aligned} F(x) &= \int \limits_0^\infty e^{s (j-\frac{1}{x})} \left[ \sum_{n=0}^N \frac{s^n}{n!} \right]^{k-j} ds \\[6pt] &= \int \limits_0^\infty e^{s (j-\frac{1}{x})} \sum_{i=0}^{N(k-j)} a_i s^i \ ds \\[6pt] &= \sum_{i=0}^{N(k-j)} a_i \int \limits_0^\infty e^{s (j-\frac{1}{x})} s^i \ ds \\[6pt] \end{aligned} \end{equation}$$

Now, assume that $j < 1/x$ so that $(j-\tfrac{1}{x}) <0$. In this case we can use the change-of-variable $m = s |j-\frac{1}{x}|$ to get $dm = -j ds$ and we have:

$$\begin{equation} \begin{aligned} F(x) &= \sum_{i=0}^{N(k-j)} a_i \int \limits_0^\infty e^{-s |j-\frac{1}{x}|} s^i \ ds \\[6pt] &= \sum_{i=0}^{N(k-j)} a_i \cdot \frac{1}{(-j) |j - \frac{1}{x}|^i} \int \limits_0^\infty e^{-s |j-\frac{1}{x}|} (s |j-\tfrac{1}{x}|)^i \ (-j) \ ds \\[6pt] &= \sum_{i=0}^{N(k-j)} a_i \cdot \frac{1}{(-j) |j - \frac{1}{x}|^i} \int \limits_0^\infty e^{-m} m^i \ dm \\[6pt] &= \sum_{i=0}^{N(k-j)} a_i \cdot \frac{i!}{(-j) |j - \frac{1}{x}|^i}. \\[6pt] \end{aligned} \end{equation}$$

This form of the function is a finite sum, where the only serious difficulty (for large $N$ or $k-j$) is to calculate the polynomial coefficients $a_0,a_1,...,a_{N(k-j)}$. Calculation of these coefficients may be cumbersome for large degree, but once they are calculated the function is then amenable to calculation.

As mentioned in another answer deriving an exact expression is definitely hard, but we can bound this sum pretty easily from above:

$$A(N,..., N)<\frac{1}{a}\sum_{n=0}^{rN}\sum_{\sum m_i=n}a^{-n}\frac{n!}{m_1!...m_r!}=\frac{1}{a}\sum_{n=0}^{rN}\Big(\frac{r}{a}\Big)^n=\frac{\Big(\frac{r}{a}\Big)^{rN+1}-1}{r-a}=\frac{\Big(\frac{k-j}{1/x-j}\Big)^{rN+1}-1}{k-1/x}$$

since as user125932 mentions, the $m_i$'s are ranging up to $N$, and not $rN$ like the multinomial distribution would demand.

Also it is obvious that due to the fact that if we only sum up to $n\leq N$ we get an exact multinomial estimate and since all the summands are positive:

$$A(N,..., N)>\frac{1}{a}\sum_{n=0}^{N}\sum_{\sum m_i=n}a^{-n}\frac{n!}{m_1!...m_r!}=\frac{1}{a}\sum_{n=0}^{N}\Big(\frac{r}{a}\Big)^n=\frac{\Big(\frac{k-j}{1/x-j}\Big)^{N+1}-1}{k-1/x}$$

which by sandwiching reduces to the correct expression of the Laplace transform of the exponential in the limit $N\rightarrow \infty$:

$$\lim_{N\to\infty}A(N,..., N)=\frac{x}{1-kx}$$

To appreciate the complexity of this sum, let us attempt an evaluation for $r=2$. The sum is close to a binomial expansion but not exactly. We can break it up as follows (here $\Delta=m+n$): $$S_2=\sum_{m,n=0}^{N}{m+n\choose m}a^{-m-n}=\Big(\sum_{\Delta=0}^{N}\sum_{m=0}^{\Delta}+\sum_{\Delta=N+1}^{2N}\sum_{m=0}^{\Delta}-\sum_{\Delta=N+1}^{2N}\sum_{m=N+1}^{\Delta}\Big){\Delta\choose m}a^{-\Delta}$$

The first two terms in the sum can be summed from the binomial theorem and the third term is a remainder.

$$S_2=\frac{\Big(\frac{2}{a}\Big)^{2N+1}-1}{\frac{2}{a}-1}-R_2$$

By changing the order of summation we find for $R_2$

$$R_2=\sum_{m=0}^{N-1}\sum_{\Delta=0}^{N-1-m}{\Delta+N+m+1\choose\Delta}a^{-\Delta}<\sum_{m=0}^{N-1}\sum_{\Delta=0}^{\infty}{\Delta+N+m+1\choose\Delta}a^{-\Delta}=\sum_{m=0}^{N-1}(\frac{a}{a-1})^{m+N+1}=\frac{a^{N+1}}{(a-1)^N}((\frac{a}{a-1})^N-1)$$

(To be continued)