Air molecules knock a little sunlight sideways, and they do it far more to short wavelengths than to long ones. Look away from the sun and you see that sideways light, which your eyes read as sky blue. Look at the sun near the horizon and you see what is left after its light crosses about 38 times more air than at noon.
A molecule is about a thousand times smaller than a wavelength of light. The light's electric field shakes the molecule's electrons, and they re-radiate in all directions with a power that grows as the fourth power of frequency. So scattering goes as 1/λ⁴. Pick two wavelengths to compare, then see what color the whole scattered spectrum makes.
How to read this: each colored bar shows how strongly air scatters that color, compared with deep red (700 nm, a light wavelength in nanometers). Blue bars are about 4 to 6 times taller. The lines show shapes only, so compare their heights with each other, not with the left axis. λ₁ and λ₂ are the two colors you're comparing (the dashed markers).
Try this: click Violet vs blue: violet scatters even more, yet the sky isn't violet. See the green curve, which weighs each color by how bright it looks to your eye.
The scattered spectrum (gold) keeps rising toward violet, so why isn't the sky violet? Three reasons. Sunlight carries less violet than blue to begin with: the sun's output peaks near 500 nm, and real sunlight dips further below about 450 nm than this smooth 5800 K curve. Ozone high in the atmosphere absorbs the ultraviolet beyond violet, so none of it joins the mix. And your eyes respond weakly at 400 nm, as the green curve shows.
The eye also never sees one wavelength. It sees the whole spread, and that blend lands on a pale blue whose dominant wavelength is near 475 nm.
import numpy as np h, c, k = 6.626e-34, 2.998e8, 1.381e-23 lam = np.arange(380, 781, 5.0) # nm, 5 nm steps def planck(l, T=5800): # sunlight as a 5800 K blackbody m = l * 1e-9 return m**-5 / np.expm1(h*c / (m*k*T)) def lobe(l, mu, s1, s2): # piecewise Gaussian return np.exp(-0.5 * ((l - mu) / np.where(l < mu, s1, s2))**2) def cmf(l): # CIE 1931 2° fit, Wyman, Sloan & Shirley 2013 x = 1.056*lobe(l,599.8,37.9,31.0) + 0.362*lobe(l,442.0,16.0,26.7) - 0.065*lobe(l,501.1,20.4,26.2) y = 0.821*lobe(l,568.8,46.9,40.5) + 0.286*lobe(l,530.9,16.3,31.1) z = 1.217*lobe(l,437.0,11.8,36.0) + 0.681*lobe(l,459.0,26.0,13.8) return np.array([x, y, z]) M = np.array([[ 3.2406, -1.5372, -0.4986], # XYZ to linear sRGB (D65) [-0.9689, 1.8758, 0.0415], [ 0.0557, -0.2040, 1.0570]]) def linear_rgb(S): # spectrum to linear sRGB, clip out-of-gamut, return np.clip(M @ (cmf(lam) @ S), 0, None) / WB # white-balance to sun def encode(rgb): # scale to brightest channel, sRGB gamma v = rgb / rgb.max() return np.where(v <= 0.0031308, 12.92*v, 1.055*v**(1/2.4) - 0.055) def dominant_nm(S, w=(0.3127, 0.3290)): # hue angle around D65 white X, Y, Z = cmf(lam) @ S a = np.arctan2(Y/(X+Y+Z) - w[1], X/(X+Y+Z) - w[0]) best = None for l in range(380, 701): x, y, z = cmf(l) d = abs(np.angle(np.exp(1j*(np.arctan2(y/(x+y+z) - w[1], x/(x+y+z) - w[0]) - a)))) if best is None or d < best[0]: best = (d, l) return best[1] sun = planck(lam) WB = np.clip(M @ (cmf(lam) @ sun), 0, None); WB /= WB.max() # sunlight shows white sky = sun * lam**-4.0 # Rayleigh: 1/λ⁴ ybar = cmf(lam)[1] print((700/450)**4) # 5.86 print(np.round(255*encode(linear_rgb(sky)))) # [129 174 255] print(dominant_nm(sky)) # 475 print(sky[lam < 450].sum() / sky.sum()) # share of scattered power print((sky*ybar)[lam < 450].sum() / (sky*ybar).sum()) # share the eye registers
The sun's direct beam loses each wavelength at its own rate as it crosses the air: T = exp(−τ·m), where τ is the vertical optical depth and m is the air mass, how many vertical atmospheres the slant path is worth. Overhead m = 1. At the horizon the path is about 38 times longer, and blue is almost entirely scattered out of the beam. Lower the sun.
How to read this: the top picture is the sky as you'd see it. Below it, on the left (top on phones), the sun's light has to cross the colored line of air to reach you; a lower sun means a longer path. On the right, how much of each color survives that trip (1 = all of it). Blue drops out first.
Try this: press Play sunset and watch the blue bars vanish.
Air mass uses the Kasten and Young (1989) formula, a fit to a standard atmosphere that includes refraction near the horizon. For comparison, a uniform shell 8.43 km thick (the scale height H = RdT/g at 288 K, as in W04 · Global circulation) gives 38.9 at the horizon. The diagram exaggerates the shell's thickness so the bend is visible.
Rayleigh optical depth uses a common sea-level fit; the exponent is 4.05 because air's refractive index also changes slightly with wavelength. Haze adds an Ångström term with exponent 1.3, and β is its optical depth at 1 µm. The sky gradient is a toy single-scattering model: flat layers, one bounce, the same scattering in every direction, no ozone. The real twilight zenith stays bluer than this because ozone absorbs orange light.
# uses lam, planck, linear_rgb, encode from Lab 1 sun = planck(lam) def air_mass(h): # Kasten & Young 1989, h = elevation in degrees return 1 / (np.sin(np.radians(h)) + 0.50572*(h + 6.07995)**-1.6364) def shell_air_mass(h, R=6.371e6, H=287*288/9.81): # uniform shell s, q = np.sin(np.radians(h)), R / H return np.sqrt((q*s)**2 + 2*q + 1) - q*s def tau(l, beta): # vertical optical depth at sea level um = l / 1000 return 0.0088*um**-4.05 + beta*um**-1.3 # Rayleigh + haze def direct_sun(h, beta=0.05): # Beer-Lambert return sun * np.exp(-tau(lam, beta) * air_mass(h)) def sky(v, h, beta=0.05): # single scattering, viewing elevation v t, ms, mv = tau(lam, beta), air_mass(h), air_mass(v) if abs(mv - ms) < 1e-6: mv = ms + 1e-6 return sun * mv * (np.exp(-t*ms) - np.exp(-t*mv)) / (mv - ms) K = 1.5 / linear_rgb(sky(90, 90)).max() # exposure: noon zenith sky def tone(rgb): # compress brightness, keep hue, sRGB gamma m = rgb.max() v = rgb * (1 - np.exp(-K*m)) / m if m > 0 else rgb return np.where(v <= 0.0031308, 12.92*v, 1.055*v**(1/2.4) - 0.055) print(air_mass(0), shell_air_mass(0)) # 37.9 38.9 T = np.exp(-tau(np.array([450, 700]), 0.05) * air_mass(0)) print(T) # kept at 450 and 700 nm print(np.round(255*encode(linear_rgb(direct_sun(3))))) # sun disc print(np.round(255*tone(linear_rgb(sky(90, 3))))) # zenith sky
The 1/λ⁴ law holds only for particles much smaller than the wavelength. Grow the particle and the wavelength dependence flattens. Once a particle is several times bigger than the wavelength (Mie scattering, and beyond it the large-particle limit) it scatters all colors about equally. The crossover drawn here is a smooth illustrative curve. A full Mie calculation would add detail to it.
How to read this: top, how much a particle favors blue light as it gets bigger. That's α: 4 means strongly, like air; 0 means not at all, like cloud droplets. Bottom, what that means color by color for the particle you picked.
Try this: click Air molecule, then Cloud droplet, and compare the two color chips below.
The curve gives α = 4 for x much less than 1 and falls toward 0 for x much greater than 1. Real Mie theory has wiggles between x ≈ 1 and 30, and for some sizes long waves even scatter a little more than short ones. This smooth curve leaves those out. Typical haze in this model has α near 1, close to the 1.3 often measured.
Clouds are droplets with radii around 10 µm, a size parameter above 100. They scatter every color alike, so they look white. A thick cloud's base looks gray because less light reaches it; the light that does arrive is still neutral.
Haze particles are a few tenths of a micron. They scatter only slightly more blue than red, so hazy skies turn milky white-gray and distant hills fade toward white.
Fine smoke from a wood fire has many particles well under the wavelength. Lit from the side against a dark background it looks bluish. The sun seen through it turns orange or red, the same reddening as Lab 2.
# uses lam, planck, linear_rgb, encode from Lab 1 sun = planck(lam) def exponent(r_um, x0=2.0): # ILLUSTRATIVE smooth curve, not Mie theory x = 2*np.pi*r_um / 0.55 # size parameter at 550 nm return 4 / (1 + (x / x0)**2) # 4 when x << 1, toward 0 when x >> 1 def scattered(r_um): # side-lit, against a dark background return sun * (lam / 550)**-exponent(r_um) def transmitted(r_um, tau550=1.0): # direct sun through the layer return sun * np.exp(-tau550 * (lam / 550)**-exponent(r_um)) for r in (0.00015, 0.05, 0.3, 10): # molecule, smoke, haze, droplet a = exponent(r) print(r, round(a, 2), round((700/450)**a, 2), np.round(255*encode(linear_rgb(scattered(r)))), np.round(255*encode(linear_rgb(transmitted(r)))))