Functional principal component analysis for sparse longitudinal data usually proceeds by first smoothing the covariance surface, and then obtaining an eigendecomposition of the associated covariance operator. Here we consider the use of penalized tensor product splines for the initial smoothing step. Drawing on a result regarding finite-rank symmetric integral operators, we derive an explicit spline representation of the estimated eigenfunctions, and show that the effect of penalization can be notably disparate for alternative approaches to tensor product smoothing. The latter phenomenon is illustrated with two data sets derived from magnetic resonance imaging of the human brain.