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!



Symmetric functions arising from continuous-time Markov chains
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](https://en.wikipedia.org/wiki/Hypoexponential_distribution). 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](https://en.wikipedia.org/wiki/Hypoexponential_distribution#Relation_to_the_phase-type_distribution), 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](https://en.wikipedia.org/wiki/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. 


