  We present a new method for constructing maps of the secondary temperature
fluctuations imprinted on the cosmic microwave background (CMB) radiation by
photons propagating through the evolving cosmic gravitational potential. Large
cosmological N-body simulations are used to calculate the complete non-linear
evolution of the peculiar gravitational potential. Tracing light rays back
through the past lightcone of a chosen observer accurately captures the
temperature perturbations generated by linear (the integrated Sachs-Wolfe or
ISW effect) and non-linear (the Rees-Sciama or RS effect) evolution. These
effects give rise to three kinds of non-linear features in the temperature
maps. (a) In overdense regions, converging flows of matter induce cold spots of
order 100 Mpc in extent which can dominate over the ISW effect at high
redshift, and are surrounded by hot rings. (b) In underdense regions, the RS
effect enhances ISW cold spots which can be surrounded by weak hot rings. (c)
Transverse motions of large lumps of matter produce characteristic dipole
features, consisting of adjacent hot and cold spots separated by a few tens of
Megaparsecs. These non-linear features are not easily detectable; they modulate
the ISW sky maps at about the 10 percent level. The RS effect causes the
angular power spectrum to deviate from linear theory at l~50 and generates
non-Gaussianity, skewing the one-point distribution function to negative
values. Cold spots of similar angular size, but much smaller amplitude than the
CMB cold spot reported by Cruz et al. are produced. Joint analysis of our maps
and the corresponding galaxy distribution may enable techniques to be developed
to detect these non-linear, non-Gaussian features. Our maps are available at
http://astro.dur.ac.uk/~cai/ISW
