Polynomial for Matrix Exponential For a square matrix $A$, the matrix exponential is defined as $$e^{At}=I/0!+At/1!+A^2t^2/2!+\cdots,$$
where $I$ is the identity matrix. There's a theorem (a corollary of the Cayley-Hamilton theorem) which states $e^{At}=a_0I+a_1At+a_2A^2t^2+\cdots +a_{n-1}A^{n-1}t^{n-1}$, where the $a_i$'s are functions of $t$ to be determined from $A$. There's another theorem (which enables us to find the $a_i's$) stating that if we define $r(\lambda)=a_0+a_1\lambda +\cdots + a_{n-1}\lambda^{n-1}$ and if $\lambda_i$ is an eigenvalue of $At$, then 


*

*$$e^{\lambda_i}=r(\lambda_i).$$ 
Further, if $\lambda_i$ is an eigenvalue of multiplicity $k>1$, then in addition we have 

*$$e^{\lambda_i}=r'(\lambda)\bigg\vert_{\lambda=\lambda_i}=r''(\lambda)\bigg\vert_{\lambda=\lambda_i}=\cdots = r^{(k-1)}(\lambda)\bigg\vert_{\lambda=\lambda_i}.$$


I am requesting a proof sketch of (2). I don't see how the derivatives come into the picture.
 A: You might be able to cook up a formal proof sketch from the material in the classic A S Householder's or F R Gantmacher's texts. Instead, here, if you were looking for a hint, I give a plausibility demonstration/illustration for the simplest n = 3 case, the simplest 3x3 diagonalizable matrix paradigm exhibiting the feature that intrigues you, the way we informally present it to physics or CS students, the primary consumers of applications. You may then generalize and formalize. 
To start with, you are evidently not utilizing the parameter t in a non-trivial way--nothing in your question would be different if you set t =1, which I will do right away.  Nevertheless, this parameter is very useful in the discussion which follows, so I will reintroduce it, and selectively set it to 0, or, in the end, 1, to comment on your expressions.
Consider the 3x3 matrix A with a single nonzero eigenvalue $\lambda_1$, and a double eigenvalue $\lambda_2$, so the Cayley-Hamilton theorem, $0=(A-\lambda_1) ( A-\lambda_2)^2$, constrains $A^3$ and higher powers to be expressible as quadratic polynomials in A, just like all of its matrix functions, including the exponential here. If there were no degeneracy, and we had 3 distinct eigenvalues, instead, the matrix exponential, like all matrix functions would be straightforward to  evaluate through Sylvester's formula. However, the  double eigenvalue displaying your feature needs Buchheim's extension. Let's first find the exponential.
Set your quadratic polynomial r(tA) to be
$$
e^{tA}= e^{\lambda_1 t } B + e^{\lambda_2 t}(C+tD) ~,  \tag 1 
$$
where B,C,D are t -independent quadratic functions of A , independent of the left hand side: the Frobenius covariants which project out eigenvectors corresponding to eigenvalues $\lambda_2$ and $\lambda_1$, determined as follows.
Evaluate (1) at t=0; operate once and twice with d/dt and evaluate at t=0 to obtain
$$ I=B+C \tag 2 
$$
$$
A=\lambda_1 B + \lambda_2 C +D   \tag 3   
$$
$$
A^2=\lambda_1^2 B + \lambda_2( \lambda_2 C + 2D) ~.  \tag 3
$$
Solve (2-4) in terms of A, to obtain 
$$
B=\frac{(A-\lambda_2)^2}{(\lambda_1-\lambda_2)^2}\\
C=I-B= \frac{(A-\lambda_2+ (\lambda_1-\lambda_2))(\lambda_1-A)}{(\lambda_1-\lambda_2)^2}\\
D= \frac{(A-\lambda_1)(A-\lambda_2)}{(\lambda_2-\lambda_1)}~.
$$
Note B is a projector onto the space of eigenvectors corresponding to $\lambda_1$, while C does the same for the two eigenvectors corresponding to $\lambda_2$; finally, D divides the characteristic polynomial above. 
Again, these Frobenius covariants are properties of A and not of the function we are going to examine, here the exponential. 
Thus,
$$
e^{tA}= e^{\lambda_1 t} \frac{(A-\lambda_2)^2}{(\lambda_1-\lambda_2)^2} + e^{\lambda_2t }  (A-\lambda_1) \frac{[(\lambda_2-A)(1+t(\lambda_1-\lambda_2))-(\lambda_1-\lambda_2)]}{(\lambda_1-\lambda_2)^2} ~,
$$
yielding, for t=1, your quadratic polynomial expression,
$$
r(\lambda)=   e^{\lambda_1 } \frac{(\lambda-\lambda_2)^2}{(\lambda_1-\lambda_2)^2} + e^{\lambda_2 }  (\lambda-\lambda_1) \frac{[(\lambda_2-\lambda)(1+\lambda_1-\lambda_2)-\lambda_1+\lambda_2]}{(\lambda_1-\lambda_2)^2} ~,
$$
with $r(\lambda_1)= e^{\lambda_1}$ and  $r(\lambda_2)= e^{\lambda_2}= e^{\lambda_2} C(\lambda_2)$.
Now, further observe directly that $r'(\lambda_2)=r(\lambda_2)$, whereas this is not the case for $r'(\lambda_1)$, that is $r'(\lambda)\neq r(\lambda.)$. This is but a demonstration, although it might well rise to a hint!



*

*PS: I realize that much of this is overkill. If you assume you have 3 distinct eigenvalues, $\lambda_1,\lambda_2, \lambda_2 +\epsilon$, and compute the  uncomplicated Sylvester formula, $r(\lambda)=   e^{\lambda_1 } \frac{(\lambda-\lambda_2)(\lambda-\lambda_2-\epsilon) }{(\lambda_1-\lambda_2)(\lambda_1-\lambda_2-\epsilon)} +  e^{\lambda_2 } \frac{(\lambda-\lambda_1)(\lambda-\lambda_2-\epsilon)}{(\lambda_2-\lambda_1)(-\epsilon)} +  e^{\lambda_2+\epsilon } \frac{(\lambda-\lambda_1)(\lambda-\lambda_2)}{(\lambda_2+\epsilon-\lambda_1) (\epsilon)}$,  obtaining the above expression to lowest (zero) order in ε; whence $r(\lambda_i)$ for all 3 i's and also $r'(\lambda_i)$, and only then taking the $\epsilon\to 0$ limit, you appreciate the inevitability of your statement 2. at $\lambda_2$.

*It is a generic feature of Lagrangian interpolation. This generalizes to higher degeneracies, mutatis mutandis. Do you really want me to carry it out more explicitly? 
