Infant cortical surface templates play an essential role in spatial normalization of cortical surfaces across individuals in pediatric neuroimaging analysis. However, existing infant surface templates have two major limitations in functional MRI analysis. First, they are constructed by co-registration of cortical surfaces based on structural attributes, which cannot lead to accurate functional alignment, due to the highly variable relationship between cortical folds and functions. Second, they are constructed by simply averaging co-registered cortical attributes, which is sensitive to registration errors and lead to blurred attribute patterns on templates, thus deteriorating the accuracy in spatial normalization. Therefore, construction of infant cortical functional templates encoding sharp functional architectures is critical for infant fMRI analysis. To this end, we construct the first set of spatiotemporal infant cortical surface functional templates using Wasserstein barycenter and a state-of-the-art functional feature, namely the gradient density of functional connectivity. To address the first issue, we leverage functional gradient density to drive surface registration to improve inter-individual functional correspondences. To address the second issue, we compute templates based on the Wasserstein barycenter of functional gradient density maps across individuals. The motivation is that Wasserstein barycenter represents a meaningful mean under the Wasserstein distance metric, which takes into account the alignment of local spatial distribution of cortical attributes and thus is robust to registration errors, leading to sharp and detailed patterns on templates. Experiments on a dataset with 207 fMRI scans between 0 and 2 years of age show the validity and accuracy of our constructed infant cortical functional templates.