  We present an in-depth and systematic variability study of a sample of 20
powerful blazars, including 12 BL Lacs and 8 flat spectrum radio quasars,
applying various analysis tools such as flux distribution, symmetry analysis,
and time series analysis on the decade-long Fermi/LAT observations. The results
show that blazars with steeper \gama-ray spectral indexes are found to be more
variable; and the \gama-ray flux distribution closely resembles lognormal
probability distribution function. The statistical variability properties of
the sources as studied by power spectral density analysis are consistent with
\emph{flicker noise} ($P(\nu)\propto1/\nu$) -- an indication of long-memory
processes at work. Statistical analysis of the distribution of flux rise and
decay rates in the light curves of the sources, aimed at distinguishing between
particle acceleration and energy dissipation timescales, counter-intuitively
suggests that both kinds of rates follow a similar distribution and the derived
mean variability timescales are in the order of a few weeks. The corresponding
emission region size is used to constrain location of \gama-ray production
sites in the sources to be a few parsecs. Additionally, using Lomb-Scargle
periodogram and weighted wavelet z-transform methods and extensive Monte Carlo
simulations, we detected year timescale quasi-periodic oscillations in the
sources S5 0716+714, Mrk 421, ON +325, PKS 1424-418 and PKS 2155-304. The
detection significance was computed taking proper account of the red-noise and
other artifacts inherent in the observations. We explain the results in the
light of current blazar models with relativistic shocks propagating down the
jet viewed close to the line of sight.
