Derivative of trace function $\DeclareMathOperator{\tr}{tr}$
Let $A,B$ be self-adjoint matrices and $f$ be a real differentiable function on $\mathbb{R}$ with derivative $f'$. Then why is it true that
$$
\left.\ \frac{d}{dt}\right|_0 \tr f(A+tB)=\tr (f'(A)B)
$$
This is used in the Klein's inequality. However, I'm not sure why exactly this is true in general. It's pretty clear why it's true for polnomials since we can use the commutation relation of the trace function, but it's harder to justify in general. I also checked the linked reference (E. Carlen, Trace Inequalities and Quantum Entropy: An Introductory Course, Contemp. Math. 529 (2010) 73–140) with no luck, as the author didn't give much explanation.
EDIT: After some further thought, let me provide an incomplete proof of what I got so far. Hopefully someone with better knowledge can finish the proof.
For simplicity, let $\lambda_i(A)$ denote the eigenvalues of $A$ in descending order, i.e., $\lambda_1(A) \ge \cdots \ge \lambda_d (A)$. Then
$$
\tr \left( \frac{f(A+tB)-f(A)}{t}\right) = \sum_i \frac{1}{t}[f(\lambda_i(A+tB)-f(\lambda_i(A))]
$$
Notice that by Weyl's inequality (stability of eigenvalues), we see that $|\lambda_i(A+tB)-\lambda_i(A)|\le t||B||$. Hence, using an $\epsilon,\delta$ arguement, we can replace the above with
$$
\sum_i \frac{1}{t}(\lambda_i(A+tB)-\lambda_i(A)) f'(\lambda_i(A))
$$
Now first assume that $A$ has a simple spectrum, then $A+tB$ is also simple for sufficiently small $t$. Then by Hadarmard's variation formula, we see that
$$
\frac{1}{t}(\lambda_i(A+tB)-\lambda_i(A)) \to \langle i|B| i\rangle
$$
where $|i\rangle$ is the corresponding eigenvector (unique up to phase since we are assuming that $A$ is simple) to $\lambda_i(A)$. Plugging all this back in, we see that the formula at least holds when $A$ is simple.
EDIT 2. I think I now have a way of dealing with degenerate eigenvalues. I will provide a sketch and fill in the details later (if someone else doesn't point out an error).
Let $\lambda_1 (A)=\cdots =\lambda_r(A)$ be the degenerate eigenvalues. Then for sufficiently small $t$, the eigenvalues $\lambda_i (A+tB),i=1,...,r$ will not touch the other eigenvalues (Weyl's inequality again). Let us use the Riesz projector
$$
P_A =\frac{1}{2\pi i} \oint_\Gamma \frac{dz}{A-z}
$$
where $\Gamma$ is some "smooth" contour around the $\lambda_1 (A)=\cdots =\lambda_r(A)$ and its interior does not contain any other eigenvalues. By Weyl's inequality, we can assume that $\lambda_i(A+tB),i=1,...,r$ are still in the interior of $\Gamma$ for sufficiently small $t$. Notice that
$$
\frac{d}{dt} \Big|_0 \tr {((A+tB)P_{A+tB})} = \tr(BP_A)
$$
where I got some inspiration from @Ruy's comment and used the fact that
\begin{align}
\frac{d}{dt}\Big|_0 \tr{(A(P_{A+tB}-P_A))}&=\tr A\oint_\Gamma \frac{dz}{(z-A)^2}B \\
&= \sum_{i=1}^r \oint_\Gamma \lambda_i(A)\frac{1}{(z-\lambda_i(A))^2} dz \langle i|B|i\rangle \\
&=0
\end{align}
Hence, if we combine this with the previous part, we see that the equality holds.
My proof is a little convoluted, so I would still hope to see a more straightforward approach
 A: Lemma 1.  Let $f$ and $g$ be real valued functions of class $C^1$ defined on a neighborhood $(a-\delta , a+\delta )$
of the real number $a$, such that
$$
  f(a)=g(a),\quad \text{and}\quad  f'(a)=g'(a).
  $$
Also let $\lambda :U\to (a-\delta , a+\delta )$ be any function defined on a neighborhood $U$  of zero
such that
$$
  |\lambda (t)-a|\leq c|t|, \quad \forall t\in  U,
  $$
where $c$ is a given  positive  constant.
Then
$$
  \lim_{t\to 0} \frac{f(\lambda _t) - g(\lambda _t)}t = 0.
  $$
Proof.  By the mean value Theorem,  for each $t$ in $U$, we may write
$$
  f(\lambda _t) = f(a) + f'(\xi _t)(\lambda _t-a),
  $$
and
$$
  g(\lambda _t) = g(a) + g'(\eta _t)(\lambda _t-a),
  $$
where $\xi _t$ and $\eta _t$ lie in between $\lambda _t$ and $a$.  Therefore
$$
  |f(\lambda _t) - g(\lambda _t)| = $$ $$ =
  |f'(\xi _t)-g'(\eta _t)||\lambda _t-a| \leq
  c|f'(\xi _t)-g'(\eta _t)| |t|.
  $$
Since both $\xi _t$ and $\eta _t$ converge to $a$,  as $t\to 0$,  we get
$$
  \lim_{t\to 0} \left|\frac{f(\lambda _t) - g(\lambda _t)}t\right| \leq $$ $$ \leq
  \lim_{t\to 0} c|f'(\xi _t)-g'(\eta _t)| =
  c|f'(a)-g'(a)| = 0.
  $$
QED
Lemma 2.  Let $A$ and $B$ be $n\times n$ self-adjoint complex matrices and
let $f$ and $g$ be real valued functions of class $C^1$ defined on a neighborhood of $\sigma (A)$, such that
$$
  f(a)=g(a),\quad \text{and}\quad  f'(a)=g'(a), \quad \forall a\in \sigma (A).
  $$
Then
$$
  \lim_{t\to 0} \frac{f(A+tB)-g(A+tB)}t = 0.
  $$
Proof.  Let $a_1\leq a_2\leq \cdots \leq a_n$ be the eigenvalues of $A$, and let $\lambda _1(t)\leq \lambda _2(t)\leq \cdots \leq \lambda _n(t)$ be the eigenvalues of $A+tB$.
By Weyl's inequality (stability of eigenvalues), we have that
$$
  |\lambda _i(t)-a_i|\leq |t| \|B\|.
  $$
So
$$
  \lim_{t\to 0}\left\|\frac{f(A+tB)-g(A+tB)}t\right\| = $$ $$ =
  \lim_{t\to 0}\max_{1\leq i\leq n}  \left|\frac{f(\lambda _i(t))-g(\lambda _i(t))}t\right| = 0,
  $$
by
Lemma (1).  QED
Theorem 3. Let $A$ and $B$ be $n\times n$ self-adjoint complex matrices and
let $f$ be a real valued function of class $C^1$ defined on a neighborhood of $\sigma (A)$.  Then
$$
  \frac d{dt}\Big|_{t=0} \text{tr}(f(A+tB)) = \hbox{tr}(f'(A)B).
  $$
Proof.  Assume first that $f$ admits a  holomorphic extension to  a complex neighborhood of the spectrum of $A$.   Then
$$
  \frac d{dt}\Big|_{t=0} f(A+tB) =
  \frac d{dt}\Big|_{t=0} \frac{1}{2\pi}\oint f(z)(z-A-tB)^{-1}\,d z =
  $$
$$
  =   \frac{1}{2\pi}\oint f(z)(z-A)^{-1} B(z-A)^{-1}\,d z
  $$
... which has the same trace as ...
$$
  \frac{1}{2\pi}\oint f(z)(z-A)^{-2} B\,d z =  $$ $$ =
 \left(\frac{1}{2\pi}\oint f(z)(z-A)^{-2} \,d z\right) B  = f'(A)B.
  $$
Back to  the general case,  let $p$ be a real polynomial such that
$$
  p(a)=f(a),\quad \text{and}\quad  p'(a)=f'(a), \quad \forall a\in \sigma (A).
  $$
We then have that
$$
  \left\|\frac{f(A+tB) - f(A)}t  -\frac{p(A+tB) - p(A)}t   \right\| = $$ $$ =
  \left\|\frac{f(A+tB) - p(A+tB)}t   \right\|,
  $$
which converges to zero by Lemma (2).  Taking the limit as $t\to 0$,  it  then follows that
$$
  \frac d{dt}\Big|_{t=0} \text{tr}(f(A+tB)) =
  \frac d{dt}\Big|_{t=0} \text{tr}(p(A+tB)) = $$ $$ =
  \hbox{tr}(p'(A)B) =
  \hbox{tr}(f'(A)B).
  $$
QED
