  We study the gamma-ray variability of 13 blazars observed with the Fermi
Large Area Telescope (LAT). These blazars have the most complete light curves
collected during the first 4 years of the Fermi sky survey. We model them with
the Ornstein-Uhlenbeck (OU) process or a mixture of the OU processes. The OU
process has power spectral density (PSD) proportional to 1/f^alpha with alpha
changing at a characteristic time scale, tau_0, from 0 (tau>>tau_0) to 2
(tau<<tau_0). The PSD of the mixed OU process has two characteristic time
scales and an additional intermediate region with 0<alpha<2. We show that the
OU model provides a good description of the Fermi/LAT light curves of three
blazars in our sample. For the first time we constrain a characteristic
gamma-ray time scale of variability in two BL Lac sources, 3C 66A and PKS
2155-304 (tau_0=25 day and tau_0=43 day, respectively, in the observer's
frame), which are longer than the soft X-ray time scales detected in blazars
and Seyfert galaxies. We find that the mixed OU process approximates the light
curves of the remaining 10 blazars better than the OU process. We derive limits
on their long and short characteristic time scales, and infer that their
Fermi/LAT PSDs resemble a power-law function. We constrain the PSD slopes for
all but one source in the sample. We find hints for sub-hour Fermi/LAT
variability in four flat spectrum radio quasars. We discuss the implications of
our results for theoretical models of blazar variability.
