We propose two methods based on the functional principal component analysis (FPCA) to estimate smooth derivatives for a sample of observed curves with a multidimensional domain. We apply the eigendecomposition to a) the dual covariance matrix of the derivatives; b) the dual covariance matrix of the observed curves, and take derivatives of their eigenfunctions. To handle noisy and discrete observations, we rely on local polynomial regression. We show that if the curves are contained in a finite-dimensional function space, the second method performs better asymptotically. We apply our methodology in simulations and an empirical study of option implied state price density surfaces. Using call data for the DAX 30 stock index between 2002 and 2011, we identify three components that are interpreted as volatility, skewness and tail factors, and we find evidence of term structure variation.