Aims: This work aims to provide a theoretical formulation of surface brightness fluctuations (SBF) in the framework of probabilistic population synthesis models that have no deterministic relations between the different stellar components of a population but only relations on average, and to distinguish between the different distributions involved in the definition of SBF. Methods: By applying the probabilistic theory of stellar population synthesis models, we estimate the shape (mean, variance, skewness, and kurtosis) of the distribution of fluctuations across resolution elements, and examine the implications for SBF determination, definition and application. Results: We distinguish between three definitions of SBF: (i) stellar population SBF, which can be computed from synthesis models and provide an intrinsic metric for fit for stellar population studies; (ii) theoretical SBF, which include the stellar population SBF plus a term accounts for the distribution of the number of stars per resolution element ψ(N); theoretical SBF that coincides with the Tonry & Schneider (1998) definition in the special case when ψ(N) has a Poisson distribution. We find that the Poisson contribution to theoretical SBF is around 0.1 stellar population SBF and is negligible; (iii) observational SBF. We present alternative ways to compute the SBF and extend the application of stellar population SBF to defining a metric for fitting for standard stellar population studies. Conclusions: We demonstrate that SBF are observational evidence of a probabilistic paradigm in population synthesis, where integrated luminosities have an intrinsic distributed nature, and they rule out the commonly assumed deterministic paradigm of stellar population modeling.