Differential calculus of functions of self-adjoint operators Let $H$ be a Hilbert space over $\mathbb{C}$. Fix a self-adjoint operator $A:D(A)\rightarrow H$ and a Borel function $f:\mathbb{R}\rightarrow\mathbb{C}$. The operator $f(A)$ is defined by the spectral theorem.
I wonder whether under certain assumptions on $f$ and $A$, something of the following form $$f(A+B）=f(A)+L(B)+o(\|B\|)$$ can be concluded for $B$ bounded and symmetric. Here $L$ denotes a bounded ($\mathbb{R}$-)linear map from the space of (symmetric) bounded operators to the space of bounded operators on $H$. (In general, self-adjointness is lost under scaling by a complex number, so I can only hope to deal with real scalars.)
In my specific problem, $H$ is finite-dimensional, $f$ is smooth but not necessarily analytic. I am also interested in what happens for infinite-dimensional $H$ (and thus possibly stronger assumptions on $f$.)
 A: J. T. Schwartz (1955) has something like this, for analytic $f$. He says a function $f(A)$, defined for $A$ in some open set of bounded operators on a Banach space, has Frechet derivative $L$ (a linear map from bounded operators to bounded operators, which he denotes $df(A)/dA$), if
$$
\left.\frac d{d\epsilon}f(A+\epsilon B)\right|_{\epsilon=0} = L(B),\qquad\text{i.e.}\qquad
f(A+\epsilon B)=f(A)+\epsilon L(B)+o(\epsilon).
$$
He proves (p. 375):

Theorem 2: Let $f(z)$ be a function analytic in a neighborhood $D$ of the spectrum $\sigma(A)$ of a bounded operator $A$. Then $f(S)$ is defined for every operator $S$ in a neighborhood $U$ of $A$, and the operator-valued function $f(S)$ of $S$ has a Frechet derivative $df(A)/dA$.

A: Assume that $\|A\|=\|B\|=1$, just to fix some normalisation. Our aim is to show that the limit $\frac{f(A+\varepsilon B) - f(A)}{\varepsilon}$ exists and is linear in $B$. By using a mollifier, we may assume that $f$ is compactly supported since we only care about values of around the spectrum, which is contained in the interval $[-2,2]$. It follows that the derivative is a Fourier transform of a nice function $g$, in particular integrable. We may then write
$$
f(x)-f(y) = i \int_{\mathbb{R}} \frac{e^{-isx}-e^{-isy}}{s} g(s) ds.
$$
Plugging in $A+\varepsilon B$ for $x$ and $A$ for $y$ we get
$$
f(A+\varepsilon B) - f(A) = i \int_{\mathbb{R}} \frac{e^{-is(A+\varepsilon B)} - e^{-isA}}{s} g(s) ds.
$$
Note now that the following formula holds (which essentially follows from the fact that it is easy to differentiate exponentials):
$$
e^{-iS} - e^{-iT} = -i\int_{0}^{1} e^{itS} (S-T) e^{i(1-t)T} dt.
$$
If you apply it for $S= s(A+\varepsilon B)$ and $T=sA$, you get
$$
f(A+\varepsilon B) - f(A) = \int_{0}^{1} \left(\int_{\mathbb{R}} e^{-its(A+\varepsilon B)} \varepsilon B e^{-i(1-t)sA}g(s)ds\right) dt,
$$
so
$$
\frac{f(A+\varepsilon B) - f(A)}{\varepsilon} = \int_{0}^{1} \left(\int_{\mathbb{R}} e^{-its(A+\varepsilon B)} B e^{-i(1-t)sA}g(s)ds\right) dt.
$$
We can now use dominated convergence theorem to show that the limit exists and is equal to
$$
\lim_{\varepsilon\to 0} \frac{f(A+\varepsilon B) - f(A)}{\varepsilon} = \int_{0}^{1} \left(\int_{\mathbb{R}} e^{-itsA} B e^{-i(1-t)sA}g(s)ds\right) dt.
$$
To conclude, one exploits Fourier transform to write our function $f$ as an integral of analytic functions for which the result is easy.
EDIT: Once we know the formula for the derivative, we can actually check that this is a Fréchet derivative. The formula for the derivative at point $A$ in the direction $B$ is
$$
D_{A}(B):=\int_{0}^{1} \left(\int_{\mathbb{R}} e^{-itsA} B e^{-i(1-t)sA}g(s)ds\right) dt.
$$
So we get 
$$
f(A+B) - f(A) - D_{A}(B) = \int_{0}^{1} \left(\int_{\mathbb{R}} (e^{-its(A+B)}-e^{-itsA}) B e^{-i(1-t)sA}g(s)ds\right) dt.
$$
From the formula $e^{-iS} - e^{-iT} = -i\int_{0}^{1} e^{itS} (S-T) e^{i(1-t)T} dt$ we conclude that $\|e^{-iS} - e^{-iT}\| \leqslant \|S-T\|$, so
$$
\|f(A+B) - f(A) - D_{A}(B)\| \leqslant \int_{0}^{1} \left(\int_{\mathbb{R}} ts \|B\|^2 |g|(s) ds \right) dt,
$$
from which we get
$$
\frac{\|f(A+B) - f(A) - D_{A}(B)\|}{\|B\|} \leqslant \|B\| \int_{0}^{1} t \left(\int_{\mathbb{R}} s |g|(s) ds\right) dt.
$$
So we only need to assume additionally that $\int_{\mathbb{R}} s |g|(s) ds < \infty$, which is satisfied under mild smoothness assumption on $f$.
A: Suppose that $A$ has discrete spectrum consisting of   eigenvalues with finite multiplicities 
$$ 0<\lambda_1 < \lambda_2<\cdots   $$
with $\lambda_n\to\infty$ as $n\to \infty$. Denote by $(\psi_n)$   an orthonormal  basis consisting of e-vectors of $A$.  In this basis $A$ is an (infinite) diagonal matrix $\DeclareMathOperator{\diag}{diag}$
$$ A=\diag(\lambda_1,\lambda_2,\dotsc, ...).$$
Then $\newcommand{\ve}{\varepsilon}$
$$f(A+\ve I)=\diag( f(\lambda_1+\ve), f(\lambda_2+\ve),\dotsc),..). $$
Assume that $\newcommand{\bR}{\mathbb{R}}$  $f:\bR\to\bR$ is $C^2$. Then
$$ f(\lambda_n+\ve I)=f(\lambda_n)+\ve f'(\lambda_n)\ve+ \frac{1}{2}f''(\xi_n)\ve^2,\;\;\xi_n\in (\lambda_n,\lambda_n+\ve). $$
Suppose more concretely that $f(x)=x^3$. The above shows that $f(A+\ve I)-f(A)$ is not bounded.
This may not be as surprising since $A^3$ is not a bounded operators.  Here is a more interesting examples. 
Suppose for simplicity that  $\lambda_n=n$ and $f$ is a  $C^2$ function such that for $|x-n|<0.1$ we have $f(x)=n^2\cos (x-n)-n^2+1$. Note that 
$$f(\lambda_n)=1,\;\;\forall n $$
so $f(A)=I$.  Consider the bounded operator
$$ B=\diag( 1^{-1/2}, 2^{-1/2},\dotsc, n^{-1/2},\dotsc, ). $$
Then
$$f(A+\ve B)-f(A)=\diag(\dotsc, f(n+\ve n^{-1/2})-f(n),\dotsc), $$
and   we observe  that,  if $\ve<0.1$, then
$$   f(n+\ve n^{-1/2})-f(n) =n^2\cos\ve n^{-1/2}.  $$
We deduce that $f(A+\ve B)-f(A)$ is not bounded.
Remark  (a) Here is a possible reformulations.   There are  several natural topologies on the space of  closed selfadjoint operators; see this paper.  One of them, called  gap topology in the above paper  is defined by a certain metric $\gamma$ (distance between the graphs)   and has the property
$$\gamma(A_n,A)\to 0 \,\Longleftrightarrow\; \Vert f(A_n)- f(A)\Vert\to 0,\;\;\forall f\in C_0(\bR), $$
where $C_0(\bR)$ denotes the space of   continuous functions $f:\bR\to\bR$ such that 
$$\lim_{t\to\pm\infty} f(t)=0. $$
It may be the  case that if $f\in C^2(\bR)$ is such that $f,f', f''\in C_0(\bR)$ then a conclusion of the type you've formulated could be true.
(b) Let me  mention a closely related question. The set of closed selfadjoint Fredholm operators on $H$, equipped with the above gap topology  can be  organized as a   Banach manifold.  More precisely it is an open dense subset $\newcommand{\eO}{\mathscr{O}}$ $\eO$ of  $\DeclareMathOperator{\Lag}{Lag}$ $\Lag(H\oplus L)$  Grassmannian of Lagrangian subspaces of $H\oplus H$; see this paper or  this paper  for more details.  The above remark  shows that any $f\in C_0(\bR)$ defines a continuous function $\hat{f}:\eO\to\bR$, $A\mapsto f(A)$ and I am $52$% sure  that it extends to $\Lag(H\oplus H)$. Is it true that if $f\in C_0(\bR)$  is such that $f'\in C_0(\bR)$ that $\hat{f}$ is a $C^1$-function on $\eO$?
