Hcfyk

Assume that $f:mathbb Rtomathbb R$ is continuous.
Given a real symmetric matrix $Aintext{Sym}(n)$, we can define $f(A)$ by applying $f$ to its spectrum. More explicitly,
$$ f(A):=sum f(lambda)P_lambda,qquad A=sumlambda P_lambda. $$
Here both sums are finite, and the second one is the decomposition of $A$ as a linear combination of orthogonal projections ($P_lambda$ is the projection onto the eigenspace for the eigenvalue $lambda$, so that $P_lambda P_{lambda'}=0$). Such decomposition exists and is unique by the spectral theorem.

I guess it is well known that $f:text{Sym}(n)totext{Sym}(n)$ is continuous.

Assuming $fin C^infty(mathbb R)$, is the induced map $f:text{Sym}(n)totext{Sym}(n)$ also smooth?

I think I can show that it is (Fréchet) differentiable everywhere, but I am wondering whether it is always $C^1$ or even $C^infty$.

Yes. The can be derived from the resolvent formalism.

I'll just do the $C^1$ case and leave higher derivatives as an exercise - ask if it's not clear how to generalize. I am basically using formula (2.7) of "On differentiability of symmetric matrix valued functions", Alexander Shapiro, http://www.optimization-online.org/DB_HTML/2002/07/499.html. (I think there's a typo in the middle case of the display preceeding (2.7): $f(mu_j)/(mu_j-mu_k)$ should be $(f(mu_j)-f(mu_k))/(mu_j-mu_k).$) Shapiro's paper references "(cf., [4])", where [4] is the 600-page textbook "Perturbation Theory for Linear Operators" by T. Kato, but I don't know if that is helpful for this specific question.

I will call the induced map $f^*$ to distinguish it from $f.$ I'll also call the dimension $p$ instead of $n.$

It suffices to show $f^*$ is $C^1$ for matrices with eigenvalues in a given bounded interval $J.$ Approximate $f$ by polynomials $f_n$ such that $sup_{xin J}|f(x)-f_n(x)|to 0$ and $sup_{xin J}|f'(x)-f_n'(x)|to 0.$ Since $f_n$ is analytic, $f^*_n$ can be evaluated using resolvents:

$$f_n^*(X) = frac{1}{2pi i}int_C f_n(z)(z I_p - X)^{-1} dz$$
where $C$ is an anticlockwise circle in the complex plane with $J$ in its interior. For $Hinmathrm{Sym}(p),$

begin{align*}
f_n^*(X+H)
&= frac{1}{2pi i}int_C f_n(z)(z I_p - X-H)^{-1} dz\
&= frac{1}{2pi i}int_C f_n(z)(z I_p - X)^{-1}+f_n(z)(z I_p - X)^{-1}H(z I_p - X)^{-1} +dots dz\
&= frac{1}{2pi i}int_C f_n(z)sum_{lambda}(z-lambda)^{-1}P_lambda +f_n(z)sum_{lambda_1,lambda_2}(z-lambda_1)^{-1}(z-lambda_2)^{-1}P_{lambda_1}HP_{lambda_2}+dots dz\
&= f_n^*(X)+sum_{lambda_1,lambda_2} P_{lambda_1} H P_{lambda_2}int_0^1 f'_n(tlambda_1+(1-t)lambda_2)+dots dt
end{align*}

The second equality uses the Taylor expansion $$(A-H)^{-1}=A^{-1}+A^{-1}HA^{-1}+dots$$ with $A=z I_p-X.$
The third equality uses $(zI_p - X)^{-1}=sum_lambda (z-lambda)^{-1} P_lambda.$ The fourth equality uses $int_C f_n(z)(z-lambda)^{-1}(z-mu)^{-1}dz =int_0^1 f'_n(tlambda+(1-t)mu)dt.$

This gives a bound

$$|Df^*_n(X)H| leq c_p|H|cdot sup_{xin J}|f'_n(x)-f_n'(x)|$$

for some constant $c_p>0,$ where $|cdot|$ is any matrix norm. This shows that $f^*$ can be approximated arbitrarily well in the $C^1$ norm, which means it's $C^1.$

Yes. To show that f(A) is n-times differentiable at A=B, simply interpolate f by a polynomial P so that P and its derivatives up to order n agree with f on the spectrum of B. Clearly P(A) is n-times differentiable, and it isn't too much work to show that f and P have the same nth derivative.