Does sky luminance really decrease in steps as the Sun goes deeper under the horizon? Playing with some atmospheric scattering simulations, I've come across a fact that, as the Sun goes lower under the horizon, sky luminance (neglecting sources of light other than the Sun) appears to not decrease very smoothly – instead it decreases in a sort of steps. See e.g. the plot of zenith luminance with Sun elevation (all plots here assume observer on the ground):

or horizon luminance (at the Sun's azimuth, so the peak here corresponds to the solar disk):

After some experiments I've concluded that each new step appears after I add a new scattering order taken into account. If I e.g. take only 3 or 2 scattering orders into account, zenith luminance appears to look like this:

and horizon luminance:

These steps correspond to the events when each order of scattering gets fully into the Earth's shadow, so that the next order becomes dominant.
Now, what I first thought is that this jumpiness is an artifact of some discretization in the simulation. Indeed, until I increased resolution in solar elevation to insane levels, I got some bogus steps not even related to any physical events. But even after increasing resolutions in various variables (Sun elevation, view zenith angle, sun-camera angle, camera altitude) and taking really many (150) orders of scattering into account these smooth steps still remain, and begin to seem physical.
The good part of this result is that it actually is practically falsifiable. Namely, it predicts fast decrease of zenith luminance till about -7.5° Sun elevation and of horizon luminance till about -15°, after which luminance decreases considerably slower until about -18° for zenith and -25° for horizon. Although the -25° seems too deep for practical measurements (due to airglow, zodiacal light etc.), the others appear to be within the astronomical (or brighter) twilight.
So my question is, has this steppy luminance change actually been measured in real Earth atmosphere?
 A: Such pronounced steps appeared to have been an artifact of the calculations in the model. Namely, integrals of radiance over all directions were done numerically, with 16 steps in $\theta$ and in $\varphi$ spherical coordinates. This led to poor resolution (in $\theta$ direction) of the horizon, where the red sky of the sunset/sunrise is a relatively thin strip, which shrinks with Sun setting deeper. After I increased number of samples, I got very different radiances:

The 64 and 96 sample curves are nearly indistinguishable here, meaning that the computation has converged to some useful precision. Here's the 96-sample version alone, to better see its smoothness:

Some very small and smooth steps are still present after the first step, but I think they'll never be seen in reality due to the noise and other, higher-intensity light effects of the Earth's atmosphere and space.
The first step doesn't get any smoother after increasing number of samples, but it isn't much surprise either: it's the border of the normal Earth's shadow coming into zenith, which should be expected to affect further evolution of zenith radiance as the Sun sets deeper.
