Differential operator acting on a Vandermonde determinant - an identity in my endeavors I've stumbled upon the following identity:
$$
\prod_{i=1}^N \left ( 1 + \frac{\partial}{\partial a_i}\right ) \left ( 
 \Delta(a) \prod_{i=1}^N a_i \right ) = \Delta(a) \int_0^\infty dt e^{-t} \prod_{i=1}^N (t+a_i) 
$$
with Vandermonde determinant $\Delta (a) = \det \left ( a_j^{i-1} \right )_{i,j=1...N}$. I find it surprisingly hard to prove, maybe someone knows this formula or could suggest any way to derive it? I've tried to utilize Schur polynomials as the LHS is the $s_{(1_N)}$ but to no avail.
EDIT: 
My attempt of a proof started from the RHS where I use the generating functional of elementary symmetric functions:
$$
\prod_{i=1}^N (t + a_i ) = \sum_{n=0}^N \sigma_n (a) t^{N-n}
$$
with elementary symmetric functions defined as
$$
\sigma_n(a) = \sum_{1\leq i_1 < i_2 < ... < i_n \leq N} a_{i_1} a_{i_2} ... a_{i_n}
$$
which gives 
$$
\text{RHS} = \Delta(a) \sum_{n=0}^N (N-n)!\sigma_n(a) 
$$
 A: I am not sure about the following being the easiest proof but this is the one I can find.
Let us start with the lemma (that can be also proven by an expansion of the determinant)
Lemma 1 Let $D\left(a_1,a_2,\ldots,a_N\right) = \det\left\{f_{i,j}\left(a_j\right)\right\}_{i,j=1}^N$. Then
\begin{equation}
\prod\limits_{i=1}^N \left(1+\frac{\partial}{\partial a_i}\right)D\left(a_1,a_2,\ldots,a_N\right) = 
\det\left\{f_{i,j}\left(a_j\right)+f'_{i,j}\left(a_{i,j}\right)\right\}_{i,j=1}^N.
\end{equation}
Corollary 1 The LHS of the initial identity can be written as $\det\left\{a_j^{i-1}\left(a_j+i\right)\right\}_{i,j=1}^N$.
Proof of Lemma 1 If one can show
\begin{equation}
 \left(1+\frac{\partial}{\partial a_N}\right)D\left(a_1,a_2,\ldots,a_N\right) = 
\det\left\{f_{i,j}\left(a_j\right)+\delta_{N,j}f'_{i,j}\left(a_{i,j}\right)\right\}_{i,j=1}^N,
\end{equation}
then the result follows from $N$ iterations of the above. Let $M = \left\{f_{i,j}\left(a_j\right)\right\}_{i,j=1}^N$ be an $N\times N$ matrix. Then
$D = \det M$ and by the Jacobi's formula one has
\begin{equation}
\left(1+\frac{\partial}{\partial a_N}\right)D\left(a_1,a_2,\ldots,a_N\right) = 
D\left(a_1,a_2,\ldots,a_N\right)\left(1+\mathrm{Tr} \, M^{-1} \frac{\partial M}{\partial a_N}\right).
\end{equation}
The derivative $\frac{\partial M}{\partial a_N}$ is the rank one matrix with the last column being only one non-zero column, and we can write $\frac{\partial M}{\partial a_N} = 
\underbrace{\left(f'_{1,N}\left(a_N\right), f'_{2,N}\left(a_N\right), \ldots,f'_{N,N}\left(a_N\right)\right)^T}_{u}\underbrace{\left(0,0,\ldots,1\right)}_{v^T}$. This leads to
\begin{eqnarray}
\left(1+\frac{\partial}{\partial a_N}\right)D\left(a_1,a_2,\ldots,a_N\right) &= &
\det\left(M\right)\left(1+\mathrm{Tr} \, M^{-1} uv^T\right)
= \det\left(M\right)\left(1+\mathrm{Tr} \,v^TM^{-1} u\right)
\\
&=&\det\left(M\right)\left(1+v^TM^{-1} u\right) = 
\det\left(M + uv^T\right),
\end{eqnarray}
where in the last equality we have used matrix determinant lemma. Definitions of vectors $u$ and $v$ yield the result now. $\qquad \blacksquare$
We are now left with the proof of
\begin{equation}
\det\left\{a_j^{i}+ia_{j}^{i-1}\right\} = 
\Delta\left(a\right)\int\limits_{0}^{\infty}\mathrm{e}^{-t}\prod\limits_{j=1}^N\left(t+a_j\right)\,\mathrm{d}\,t.\tag{1}\label{eq:main}
\end{equation}
As usual, let $S_N$ be a set of all permutations of indexes $1,2,\ldots,N$. 
Definition For any multi-index $\mathbf{i} = \left(i_1,i_2,\ldots,i_k\right)$ with 
$1\leq i_1<i_2<\ldots<i_k\leq N$ we define a map $\tau_{\mathbf{i}}:S_N\to S_{N,k}$ by
$$
\left(\tau_{\mathbf{i}}\sigma\right)_j:=\sigma_j+\mathbf{1}_{j\in \mathbf{i}}.
$$
For an empty multi-index we define $\tau_{\emptyset} = \mathrm{id}$. For any $0\leq k\leq N$ let $S_{N,k}$ be the set of all multi-indexes $\tau_{\mathbf{i}}\sigma$ for all 
$\sigma\in S_N$ and $\mathbf{i}$ of length $k$.
Couple of obvious observations about $S_{N,k}$: this is a set of multi-indexes of length $N$, such that $\sigma = \left(\sigma_1,\sigma_2,\ldots,\sigma_N\right)$ and


*

*all $\sigma_j$ are positive integers;

*for any $1\leq i\leq N$ one has $1\leq \sigma_i \leq N+1$;

*$\sum\limits_{j=1}^N \sigma_j = \frac{N\left(N+1\right)}{2}+k$;

*it is impossible to have different $i,j,k$ such that $\sigma_i=\sigma_j=\sigma_k$.


The LHS of \ref{eq:main} can be written as
\begin{equation}
\mathrm{LHS} = \sum\limits_{\sigma \in S_N}\left(-1\right)^{\sigma}
\prod\limits_{j=1}^N a_j^{\sigma_j-1}\left(a_j+\sigma_j\right)
=\sum\limits_{\sigma \in S_N}\left(-1\right)^{\sigma}
\sum\limits_{\mathbf{i}}
\prod\limits_{j=1}^N a_j^{\left(\tau_{\mathbf{i}}\sigma\right)_j-1}
\prod\limits_{i \notin \mathbf{i}} \sigma_i.
\end{equation}
The main point of the following is a major cancellation of terms.
Proof The sum in the RHS of the above can be written as a sum over 
$\mu=\tau\sigma \in S_{N,k}$
for different $k$'s. We show that only $\mu \in S_{N,k}$ that survive are the permutations of $1,2,\ldots,N-k,N-k+2,N-k+3,\ldots,N+1$. Let $\mu \in S_{N,k}$ has all $\mu_j$ being different, than it is of the above form and $\mu=\tau_{\mathbf{i}}\sigma$ for 
$\mathbf{i} = \left(\mu^{-1}\left(N-k+2\right),\mu^{-1}\left(N-k+3\right),\ldots,\mu^{-1}\left(N+1\right)\right)$ and $\sigma = \tau^{-1}_{\mathbf{i}}\mu$. The subset of such a multi-indexes we denote $\hat{S}_{N,k}$. The coefficient of the corresponding term is equal to $\prod\limits_{i\notin \mathbf{i}} \sigma_i = \prod\limits_{i\notin \mathbf{i}} \mu_i = (N-k)!$. Now imagine $\mu$ has a pair of equal elements $\mu_p=\mu_q=r$. Then for any $\mathbf{i}, \sigma$ such that $\mu=\tau_{\mathbf{i}}\sigma$
there is a permutation $\sigma' = 
\left(
\begin{array}{cccccccc}
1 & 2 & \ldots & p & \ldots & q & \ldots & N \\
1 & 2 & \ldots & q & \ldots & p & \ldots & N 
\end{array}
\right)\sigma$ and $\mathbf{i}'$ different from $\mathbf{i}$ by changing $p$ to $q$ or otherwise (depending whether $p\in\mathbf{i}$ or $q\in\mathbf{i}$) such that $\mu=
\tau_{\mathbf{i}'}\sigma'$. Both these terms have the same coefficient, but different signs due to $\sigma$ and $\sigma'$ being different by one transposition. Therefore they cancel each other. $\quad \blacksquare$
The LHS of \ref{eq:main} can be written now as
\begin{equation}
\mathrm{LHS} = \sum\limits_{\mu\in S_{N,k}} \left(-1\right)^{\sigma}\prod\limits_{j=1}^{N}a_j^{\mu_j-1} (N-k)!.
\end{equation}
Now the RHS can be written as
\begin{equation}
\mathrm{RHS} = \left(\sum\limits_{\sigma\in S_N}
\left(-1\right)^{\sigma}\prod\limits_{j=1}^{N}a_j^{\sigma_j-1}\right)\cdot
\left(\sum\limits_{\mathbf{i}}\left(N-\left|\mathbf{i}\right|\right)!\prod\limits_{j\in\mathbf{i}} a_j\right).
\end{equation}
Using the same argument as above one can show that the only multi-indexes that will survive after multiplication have the same form as above and corresponding coefficient coincides. This proves the statement.
