Long period variable stars arise in the final stages of the asymptotic giant branch phase of stellar evolution. They have periods of up to ∼ 1000d and amplitudes that can exceed a factor of three in the I-band flux. These stars pulsate predominantly in their fundamental mode 1-3 , which is a function of mass and radius, and so the pulsation periods are sensitive to the age of the underlying stellar population 4 . The overall number of long period variables in a population is directly related to their lifetime, which is difficult to predict from first principles because of uncertainties associated with stellar mass-loss and convective mixing. The time variability of these stars has not been previously taken into account when modeling the spectral energy distributions of galaxies. Here we construct time-dependent stellar population models that include the effects of long period variable stars, and report the ubiquitous detection of this expected 'pixel shimmer' in the massive metal-rich galaxy M87. The pixel light curves display a variety of behaviors, including linearly rising and falling curves, semi-periodic curves, and sudden increases or decreases in the flux level. The observed variation of 0.1 − 1% is very well matched to the predictions of our models. The data provide a strong and novel constraint on the properties of variable stars in an old and metal-rich stellar population, and we infer that the lifetime of long period variables in M87 is shorter by approximately 30% compared to predictions from the latest stellar evolution models.In typical massive galaxies with ∼ 10 11 stars, the variation in the total light due to long period variables will be small as the summed light curves of many such stars effectively cancel each other out (with random phases the net effect scales as N −1/2 ). If the light is spread out over many (e.g., ∼ 10 4 −10 6 ) pixels, then the number of stars per pixel can range from 10 4 − 10 7 and in this regime the number of asymptotic giant branch stars per pixel is small and governed by Poisson statistics. This 'semi-resolved' regime is well known 5 and the expected surface brightness fluctuations due to Poisson statistics of rare luminous stars have been observed and studied in several hundred nearby galaxies 6-8 . We expect in this regime to be able to detect the presence of variable stars through the time dependence of the pixel flux (i.e., the pixel light curve): essentially every pixel is expected to 'shimmer' on timescales of several hundred days.In order to quantify the expected pixel shimmer, we created a stellar population model at solar metallicity that includes the time-dependent effect of long period variables. We started with a new library of stellar isochrones (Choi et al., in prep) that densely samples fast phases of stellar evolution, and assigned periods to evolved stars assuming that they pulsate in the fundamental mode 4 . We then used observations of variable stars in the Galactic bulge from the OGLE survey to estimate a period-amplitude relation in the I−band 9...