Symmetric functions arising from continuous-time Markov chains I'm interested in the following symmetric functions $s_k: \mathbb{R}_+^{k}\mapsto\mathbb{R_+}$:
\begin{align*}
s_k(x_1,x_2,\dots, x_k)&=
\int\limits_{0=t_0<t_1<\dots<t_{k-1}<t_k=1} e^{-(t_1-t_0) x_1} 
e^{-(t_2-t_1)x_2}\dots e^{-(t_k-t_{k-1})x_k}dt_1 \dots dt_{k-1} \\
&=\sum_{i=1}^k \frac{e^{-x_i}}
{\displaystyle\prod_{\substack{j=1,\dots,n\\ j\ne i}}(x_j-x_i)}
\end{align*}
The formula on the second line holds whenever the $x_i$ are all distinct; 
however it becomes degenerate if $x_i=x_j$ for some $i\ne j$. Question: is there
some "nice" pattern for the closed form expressions one obtains when some of the $x_i$ coincide? (Examples below.)
[
One motivation: consider a continuous-time Markov chain, with state space $\{1,2,3,\dots\}$, with rate $q_{ij}$ of jumping from state $i$ to state $j$ for $i\ne j$, and total rate $\lambda_i=\sum_{j\ne i} q_{ij}$ of leaving state $i$. 
Suppose the chain starts in state $1$ at time $0$. Consider the event that during the time-interval $[0,1]$, the chain follows precisely the path $1\to2\to\dots\to (k-1)\to k$ and then stays in state $k$ until the end of the interval (making $k-1$ jumps in all). 
This event has probability
$q_{12}\dots q_{(k-1)k}s_k(\lambda_1, \dots, \lambda_k)$.
]
Anyway: everything is nice and smooth and closed-form expressions for cases with repeated entries can straightforwardly be derived (using, say, L'Hôpital or whatever) - however the details appear to get messy quite quickly. I'm wondering if there is some convenient or systematic way to write those expressions (other than just calculating each case).
For example $s_2(x, y)=(e^{-x}-e^{-y})/(y-x)$,
and then $s_2(x, x)=e^{-x}$.
Taking e.g. $k=4$, we get (I think)
\begin{align*}
s_4(x,x,y,y)&=\frac{2(e^{-x}-e^{-y})}{(x-y)^3}+\frac{(e^{-x}+e^{-y})}{(x-y)^2}
\\
s_4(x,x,x,y)&=\frac{e^{-y}-e^{-x}}{(x-y)^3}-\frac{e^{-y}}{(x-y)^2}+
\frac{e^{-y}}{2(x-y)}
\\
s_4(x,y,z,z)&=\frac{e^{-x}}{(z-x)^2(y-x)} + \frac{e^{-y}}{(z-y)^2(x-y)}
\\
&\,\,\,\,\,\,\,+e^{-z}\left\{\frac{(z-x)+(z-y)+(z-x)(z-y)+z^2}{(z-x)^2(z-y)^2}\right\}.
\end{align*}
I'm aware that these things can be written in terms of multiple integrals with gamma densities; for example, if $k=k_1+\dots+k_m$ and there are $k_i$ arguments $x_i$ for each $i$, then I think
\begin{equation*}
s_k(\dots)=\int\limits_{0=t_0<t_1<\dots<t_{m-1}<t_m=1}
\prod_{i=1}^m \left\{\frac{(t_i-t_{i-1})^{k_i-1}}{(k_i-1)!} 
e^{-(t_i-t_{i-1})x_i}\right\}dt_1\dots dt_{m-1},
\end{equation*}
but I'm not sure that leads anywhere particularly nice. 
Perhaps (beguiled by the beautiful form in the case of all $x_i$ distinct)
I'm looking for something unreasonable, and these things just are what they are. Then again maybe someone has something
nice to point out!
 A: Let the random variable $t_{k-1}$ be the sum of the holding times in states $1$ up to state $k-1$.  Note that $t_{k-1}$ is a hypoexponential random variable.  Let $f_{k-1}$ denote the PDF of this random variable. From the probabilistic interpretation given by the OP, the quantity of interest can be written as:
\begin{align*}
s_k &= \frac{1}{\lambda_1 \cdots \lambda_{k-1}} \int_0^1 f_{k-1}(z) \int_{1-z}^{\infty} \lambda_k e^{-\lambda_k s} ds dz \\
&=\frac{e^{-\lambda_k}}{\lambda_1 \cdots \lambda_{k-1}} \int_0^1 f_{k-1}(z) e^{z \lambda_k} dz
\end{align*} To be sure, this basically comes from the probability that $t_{k-1}<1$ and that the holding time in state $k$ is at least $1-t_{k-1}$.
To avoid assuming that the jump rates are distinct, we will use the following matrix representation of the PDF $f_{k-1}$:
$$
f_{k-1}(z) = - \boldsymbol{\alpha}_{k-1} e^{z \Theta_{k-1}} \Theta_{k-1} \mathbf{1}_{k-1}
$$
where I borrow the notation from here, except I have introduced subscripts to $\boldsymbol{\alpha}$, $\Theta$ and $\mathbf{1}$ in order to express the number of jump rates in the hypoexponential distribution.  The matrix $\Theta_{k-1}$ is upper bidiagonal, and as long as all of the rates are positive, this matrix is invertible.  Basically, we used the fact that a hypoexponential distribution is an example of a phase-type distribution.
Let $\mathbf{I}_{k-1}$ be the $(k-1) \times (k-1)$ identity matrix.
By integrating by parts, it is straightforward to show that:
$$
\int_0^1 e^{z \Theta_{k-1}} e^{z \lambda_k} dz = (\Theta_{k-1} + \lambda_k \mathbf{I}_{k-1})^{-1} (e^{\Theta_{k-1}} e^{\lambda_k} - \mathbf{I}_{k-1} ) 
$$
Hence,
$$
s_k = \frac{-e^{-\lambda_k}}{\lambda_1 \cdots \lambda_{k-1}} \boldsymbol{\alpha}_{k-1}(\Theta_{k-1} + \lambda_k \mathbf{I}_{k-1})^{-1} (e^{\Theta_{k-1}} e^{\lambda_k} - \mathbf{I}_{k-1} ) \Theta_{k-1} \mathbf{1}_{k-1}
$$
Note that nowhere in this derivation did we assume that the jump rates are distinct.  This formula is straightforward to implement numerically.  
A: $ \text{Dear James,} $ you can write your function as a(n infinite) sum of homogeneous symmetric functions
\begin{equation*}%$
h_\ell(x_1, \dots, x_k) := \sum_{1 \leqslant i_1 \leqslant \cdots \leqslant i_\ell \leqslant k} x_{i_1} \cdots x_{i_\ell} = [t^\ell]\prod_{j = 1}^k \frac{1}{1 - t x_j}
\end{equation*}
where $ [t^\ell]f(t) $ is the $\ell$-th Fourier coefficient of $f$ (for this last equality, see wikipedia or the book by Macdonald). Now, decomposing the rational fraction, we get
\begin{equation*}%$
\prod_{j = 1}^k \frac{1}{1 - t x_j} = \sum_{j = 1}^k \frac{A_j(X)}{1 - t x_j}, \quad A_j(X) = \prod_{i \neq j} \frac{1}{1  - x_j^{-1} x_i}
\end{equation*}
thus, taking the coefficient of $ [t^\ell] $, one gets 
\begin{equation*}%$
h_\ell(x_1, \dots, x_k) = \sum_{j = 1}^k A_j(X) x_j^\ell = \sum_{j = 1}^k   \frac{x_j^{\ell- k + 1} }{ \prod_{i \neq j} (x_j - x_i) } 
\end{equation*}
Your function thus writes
\begin{equation*}%$
s_k(x_1, \dots, x_k) = \sum_{\ell \geqslant 0} \frac{(-1)^\ell}{\ell !} h_{\ell + k - 1}(x_1, \dots, x_k)
\end{equation*}
and is hence projective due to the projectivity of the $ h_i $'s. As a result, setting two variables equal amounts to get (for instance) $ h_\ell(x_1, \dots, x_{k - 1}, x_{k - 1}) $, which does not have such a simple expression. But we can continue the analysis, using the blog of Tao (https://terrytao.wordpress.com/2017/08/06/schur-convexity-and-positive-definiteness-of-the-even-degree-complete-homogeneous-symmetric-polynomials/). One can get a really nice representation of the complete homogeneous symmetric functions writing $ \frac{1}{1 - t x} = \mathbb{E}  \left( e^{ tx \mathbb{e} } \right) $ where $ \mathbb{e} \sim \mathrm{Exp}(1) $. Considering a sequence of independent such exponential random variables $ (\mathbb{e}_i)_{1 \leqslant i \leqslant k} $, one thus gets
\begin{equation*}%$
h_\ell(x_1, \dots, x_k) = [t^\ell]\prod_{j = 1}^k \frac{1}{1 - t x_j} = [t^\ell]\mathbb{E}\left( e^{ t \sum_{i = 1}^k x_i \mathbb{e}_i }  \right) = \frac{1}{\ell !} \mathbb{E}\left( \left(  \sum_{i = 1}^k x_i \mathbb{e}_i    \right)^\ell \right)
\end{equation*}
Since $ x_i \in \mathbb{R}_+ $, your function admits a representation as the Bessel expectation
\begin{equation*}%$
s_k(x_1, \dots, x_k) = \mathbb{E}\left( \sum_{\ell \geqslant 0} \frac{(-1)^\ell}{  \ell ! (\ell + k  - 1)! } \left( \sum_{i = 1}^k x_i \mathbb{e}_i    \right)^{\ell + k - 1}  \right) = \mathbb{E}\left( J_{k + 1} \left(2 \sqrt{ \sum_{i = 1}^k x_i \mathbb{e}_i  } \right)\right)
\end{equation*}
You can transform again this form using the generating series of the Bessel functions and write $ J_k(x) $ as the Fourier coefficient $ [t^k] e^{ x (t - t^{-1})/2 } $ or as a Laplace transform. Last, you can even use a stable $1/2$ random variable to write the final result as $ \mathbb{E}\left( \exp\left( Z \sum_{i = 1}^k x_i \mathbb{e}_i    \right)\right) $ where $ Z $ is a random variable independent of the $ \mathbb{e}_i $'s.
Under this form, the operation of setting two variables equal amounts to change an exponential random variable by a sum of two exponential random variables.
Best, Y.
