Modelling a large bundle of curves arises in a broad spectrum of real applications. However, existing literature relies primarily on the critical assumption of independent curve observations. In this paper, we provide a general theory for large-scale Gaussian curve time series, where the temporal and cross-sectional dependence across multiple curve observations exist and the number of functional variables, $p,$ may be large relative to the number of observations, $n.$ We propose a novel functional stability measure for multivariate stationary processes based on their spectral properties and use it to establish some useful concentration bounds on the sample covariance matrix function. These concentration bounds serve as a fundamental tool for further theoretical analysis, in particular, for deriving nonasymptotic upper bounds on the errors of the regularized estimates in high dimensional settings. As {\it functional principle component analysis} (FPCA) is one of the key techniques to handle functional data, we also investigate the concentration properties of the relevant estimated terms under a FPCA framework. To illustrate with an important application, we consider {\it vector functional autoregressive models} and develop a regularization approach to estimate autoregressive coefficient functions under the sparsity constraint. Using our derived nonasymptotic results, we investigate the theoretical properties of the regularized estimate in a "large $p,$ small $n$" regime. The finite sample performance of the proposed method is examined through simulation studies.