We model the transport of light through the atmosphere including polarization and change of wavelength as an iterated sequence of collisions, changes of wavelength, directional scatterings, and changes of the state of polarization. The collision probabilities are determined by the extinction coefficient, the change of wavelength by the relation of the collision intensities of the different `components' of the scattering medium, and the directional scattering distributions by the polarized phase function which is obtained from scattering intensities. These scattering intensities depend on the ensembles of scattering particles, the wavelength, and the polarization. We obtain them from the Mueller matrices determining the single scattering event. For easier simulation of the directional scattering the polarized phase function, which is a probability density on the unit sphere of the 3D Euclidean space, is disintegrated into the conditional polarized phase function (characterizing the azimuthal direction of scattering) and the (ordinary) phase function (characterizing the zenithal off axis direction of scattering). We characterize the path of a photon of this corpuscular stochastic process of multiple scattering by a sequence of tupels each describing the point of collision, the wavelength after collision, the direction of scattering, and the state of polarization after scattering. Starting from this model we derive exact multiple scattering lidar equations including polarization and change of wavelength. Step by step we derive from these lidar equations more specific ones: first an exact multiple scattering lidar equation with polarization (and without change of wavelength), then a multiple scattering lidar equation without polarization for collections of particles with a scattering distribution which is rotational invariant with respect to the incident beam.