Why is the sky blue?

Simple Answer

The molecules in the air scatter the light coming from the sun. However, they do not scatter the light with different frequencies equally and blue light (high frequency) gets far more scattered in the sky, giving it its color.

Coulomb's Law

To understand how this works, we are first going to need to do some refresh. This image represents the Coulomb's Law defined by:

\vec{F} = k \frac{q_1 q_2}{r^2} \vec{e_r}

Here, $\vec{F}$ represents the force applied by $q_1$ on $q_2$ , $k$ is the Coulomb's constant equal to $\frac{1}{4 \pi \varepsilon_0}$ and $\vec{e_r}$ is the unit vector going from $q_1$ to $q_2$ .

Note that if both charges are of the same sign, the force is going to be repulsive, whereas if the charges have different sign, their product will be negative and the force will be atractive, which is aligned with what we would expect.

From this force, we can derive the electric field created by $q_1$ :

\begin{aligned} \vec{F} & = q_2 \vec{E}\\ \implies \vec{E} & = k \frac{q_1}{r^2} \vec{e_r} \end{aligned}

Link Between Electric Field And Electric Potential

We also need to recall the link between the electric field and the potential in the frame of electrostatics:

\begin{aligned} V & = - \int \vec{E} . \vec{\delta l} \\ \vec{E} & = - \vec{\nabla} V \end{aligned}

Note that the electric potential $V$ is defined with an invariance to shifting by a constant, (that is why we mostly care about differences of potentials) most of the time we take that constant so that $\lim_{r \to \infty} V = 0$ . In the case of a single charged particle as before we get:

$V = k \frac{q_1}{r}$

Electrostatic Dipole

Dipole

Now we are going to see what happens when we have two opposites charges affecting the field. Note that here, $O$ is a point with no charge in the middle of $N$ and $P$ , we define $r$ as the distance between $O$ and the viewer in $M$ . Given the symmetry of the problem, the electric field has no part outside of the plane and does not depend on what plane containing the charges we are using. Thus, we can only define a single angle $\theta = (\vec{OP}, \vec{OM})$

Note that this graph is not at scale and actually represents a situation where $r >> d$ , which makes sense when we compare the distance between you and an air molecule in the atmosphere (which is going to act as the dipole) with the size of that molecule.

The electric field has the linearity propriety that makes it possible to calculate it as the sum of the two electric fields that would exist for each of both charged particle. To calculate it we are first going to compute the potential (also linear because integration is linear) and then derive the electrical field from that potential by taking the opposite of its gradient:

\begin{aligned} V & = V_+ + V_- \\ V & = kq \left( \frac{1}{PM} - \frac{1}{NM} \right) \end{aligned}

We are now going to use the fact that $r >> d$ to approximate $PM$ and $NM$ .

Dipole Approximation

Since we are very far from the dipole we are looking at, it is as we have parallel lines between the two poles and we can make the following approximation:

\begin{aligned} PM \approx P'M & = OM - OP' = r - d \cos \theta \\ NM \approx N'M & = OM + ON' = r + d \cos \theta \end{aligned}

Here $P'$ and $N'$ are, respectively, the projection of $P$ and $N$ on $(OM)$ . We can now have the potential in $M$ with the dipolar moment $p = 2qd$ :

\begin{aligned} V & = kq \left( \frac{NM - PM}{PM \cdot NM} \right) \\ & \approx kq \left( \frac{2d \cos \theta}{(r - d \cos \theta)(r + d \cos \theta)} \right) \\ & = kq \left( \frac{2d \cos \theta}{r^2 - d^2 (\cos \theta)^2} \right) \\ & \approx \frac{kp \cos \theta}{r^2} \end{aligned}

Okay, now let's compute the electric field where $\vec{e_\theta}$ is $\vec{e_r}$ rotated by $90°$ in the counter-clockwise direction:

\begin{aligned} \vec{E} & = - \vec{\nabla} V \\ & \approx - \vec{\nabla} \frac{kp \cos \theta}{r^2} \\ & = - \frac{\partial V}{\partial r} \vec{e_r} - \frac{1}{r} \frac{\partial V}{\partial \theta} \vec{e_\theta} \\ & = \frac{2kp \cos \theta}{r^3} \vec{e_r} + \frac{kp \sin \theta}{r^3} \vec{e_\theta} \\ & = \frac{kp}{r^3} (2 \cos \theta \vec{e_r} + \sin \theta \vec{e_\theta}) \end{aligned}

Rayleigh Scattering

Rayleigh scattering is, in our case, the scattering of light coming from the Sun by air molecules with a small size compared to the wavelength of the light they are scattering, i.e. $d << \lambda$ . We are going to see later why we need this hypothesis but let's see first if it is a reasonable assumption to make. Well, the visible light has wavelengths in $[400\:nm, 800\:nm]$ and the air is mainly homogeneous and composed of dinitrogen, dioxygen and water vapor which have sizes around $0.3\:nm$ . Thus, it is pretty safe to say that $d << \lambda$ .

Excitation

When the electromagnatic waves that are the light from the Sun hit the air molecules, they excite them at their frequencies and the dipoles start oscillating at those same frequencies radiating their own electromagnetic field coming from their dipolar moment.

In this case, it is as we replaced our charges $q$ by $q_0 \cos \omega t$ and $-q$ by $-q_0 \cos \omega t$ . We have the charges oscillating at the frequency of the light $f = \frac{\omega}{2 \pi}$ and here, having $d << \lambda$ allows us to be sure that the excitation is pretty homogeneous at the scale of the dipole and that the speed of the charges $v << c$ meaning that we don't have to take relativistic effects into account.

In $M$ , we have then the following retarded potential (taking into account the time we need to perceive the charges in $N$ and $P$ ):

V = kq_0 \left( \frac{\cos \omega (t - \frac{PM}{c})}{PM} - \frac{\cos \omega (t - \frac{NM}{c})}{NM} \right)

Now, using the approximation we previously made on $PM$ and $NM$ , we can simplify the cosines:

\begin{aligned} & PM = r - d \cos \theta \\ \implies & \cos \omega \left(t - \frac{PM}{c}\right) \approx \cos \left( \omega \left(t - \frac{r}{c}\right) + \frac{\omega d}{c} \cos \theta \right) \\ & = \cos \left( \omega \left( t - \frac{r}{c} \right) \right) \cos \left( \frac{\omega d}{c} \cos \theta \right) \\ & - \sin \left( \omega \left( t - \frac{r}{c} \right) \right) \sin \left( \frac{\omega d}{c} \cos \theta \right) \\ & \frac{\omega}{c} = \frac{2\pi}{\lambda} \\ \implies & \frac{\omega d}{c} = \frac{2\pi d}{\lambda} << 1 \\ \implies & \cos \omega \left(t - \frac{PM}{c}\right) \approx \cos \left( \omega \left( t - \frac{r}{c} \right) \right) - \frac{\omega d}{c} \cos \theta \sin \left( \omega \left( t - \frac{r}{c} \right) \right) \end{aligned}

In a similar way, we obtain:

\cos \omega \left(t - \frac{NM}{c}\right) \approx \cos \left( \omega \left( t - \frac{r}{c} \right) \right) + \frac{\omega d}{c} \cos \theta \sin \left( \omega \left( t - \frac{r}{c} \right) \right)

By plugging those approximations in the formula for V, we get a two kilometer long equation, which I'm going to spare you the details of, that simplifies to:

\begin{aligned} V \approx & \frac{kp_o \cos \theta}{r} \left( \frac{1}{r} \cos \omega \left( t - \frac{r}{c} \right) - \frac{\omega}{c} \sin \omega \left( t - \frac{r}{c} \right) \right) \\ \frac{\omega}{c} & = \frac{2\pi}{\lambda} >> \frac{1}{r} \\ \implies V \approx & - \frac{kp_o \omega \cos \theta}{rc} \sin \omega \left( t - \frac{r}{c} \right) \end{aligned}

It's nice, however, in electrodynamics, $V$ is not sufficient in order to derive $\vec{E}$ . We also need to compute $\vec{A}$ , the magnetic potential. Note that $\vec{e_z}$ is the unit vector in the direction $\vec{NP}$ . I will also somewhat speed up the derivation here:

\begin{aligned} \vec{I}(t) & = - q_0 \omega \sin ( \omega t ) \vec{e_{z}} \\ \implies \vec{A} & = - \frac{\mu_0 q_0 \omega}{4\pi} \int_{-d}^d \frac{\sin \omega \left( t - \frac{r(z)}{c} \right) \vec{e_z}}{r(z)} dz \\ & \approx - \frac{\mu_0 p_0 \omega}{4\pi r} \sin \omega \left( t - \frac{r}{c} \right) \vec{e_z} \\ & = - \frac{k p_0 \omega}{rc^2} \sin \omega \left( t - \frac{r}{c} \right) \vec{e_z} \end{aligned}

Electromagnetic radiation

Now that we have the electric and magnetic potentials $V$ and $\vec{A}$ , we can compute the electromagnetic field $\vec{E}$ and $\vec{B}$ radiated by the dipole:

\begin{aligned} \vec{E} & = - \vec{\nabla} V - \frac{\partial \vec{A}}{\partial t} \\ \vec{\nabla} V & = \frac{\partial V}{\partial r} \vec{e_r} + \frac{1}{r} \frac{\partial V}{\partial \theta} \vec{\theta} \\ \end{aligned}

The computation of this gradient gives two terms in $\frac{1}{r}$ and another one along $\vec{e_r}$ in $\frac{\omega}{c} = \frac{2\pi}{\lambda} >> \frac{1}{r}$ . Thus, we can approximate again:

\begin{aligned} \vec{\nabla} V & \approx \frac{kp_0\omega^2 \cos \theta}{rc^2} \cos \omega \left( t - \frac{r}{c} \right) \vec{e_r} \\ \frac{\partial \vec{A}}{\partial t} & = - \frac{kp_0\omega^2}{rc^2} \cos \omega \left( t - \frac{r}{c} \right) \vec{e_z} \\ \text{Thus, } \vec{E} & = - \frac{kp_0\omega^2}{rc^2} \cos \omega \left( t - \frac{r}{c} \right) \left( \cos \theta \vec{e_r} - \vec{e_z} \right) \\ \vec{e_z} & = \cos \theta \vec{e_r} - \sin \theta \vec{e_\theta} \\ \implies \vec{E} & = - \frac{kp_0\omega^2 \sin \theta}{rc^2} \cos \omega \left( t - \frac{r}{c} \right) \vec{e_\theta} \\ & = - \frac{\mu_0 p_0\omega^2 \sin \theta}{4\pi r} \cos \omega \left( t - \frac{r}{c} \right) \vec{e_\theta} \end{aligned}

Now we can calculate $\vec{B}$ as $\vec{\nabla} \times \vec{A}$ . Since $\vec{A}$ doesn't depend on $\varphi$ and has no component along $\vec{e_\varphi}$ , the computation of the curl in spherical coordinates is simplified to the following:

\begin{aligned} \vec{B} = \vec{\nabla} \times \vec{A} & = \frac{1}{r} \left( \frac{\partial}{\partial r}(rA_\theta) - \frac{\partial A_r}{\partial \theta} \right) \vec{e_\varphi} \\ & \approx - \frac{\mu_0 p_0\omega^2 \sin \theta}{4\pi rc} \cos \omega \left( t - \frac{r}{c} \right) \vec{e_\varphi} \end{aligned}

Here we also made another approximation of the same type as before.

Poynting Vector

Now that we have the full electromagnetic field radiated by the dipole, we can compute the Poynting vector $\vec{\Pi}$ which defines a flux density of the energy radiated by an electromagnetic field.

\begin{aligned} \vec{\Pi} & = \frac{\vec{E} \times \vec{B}}{\mu_0} \\ & = \frac{\mu_0 p_0^2 \omega^4 \sin^2 \theta}{16 \pi^2 r^2 c} \cos^2 \omega \left( t - \frac{r}{c} \right) \vec{e_r} \end{aligned}

Which averaged over time gives (since $<\cos^2> = \frac{1}{2}$ ):

\begin{aligned} <\vec{\Pi}> & = \frac{\mu_0 p_0^2 \omega^4 \sin^2 \theta}{32 \pi^2 r^2 c} \vec{e_r} \end{aligned}

Notice that the radiated power density varies in $\sin^2\theta$ , the power is zero along the axis of the dipole and is maximal along the equator. We can see on the following graph the radiated power as the radius for a given angle where the dipole is oscillating along the top-bottom axis

Dipole Radiation Graph

So, in order to get the power radiated, we need to integrate the density over a sphere of radius $r$ :

\begin{aligned} P & = \int_{S(r)} <\vec{\Pi}> . \vec{dS} \\ & = \frac{\mu_0 p_0^2 \omega^4}{32 \pi^2 c} \int_{S(r)} \left[ \frac{\sin^2 \theta}{r^2} (r d\theta) (r \sin (\theta) d\varphi) \right] \\ & = \frac{\mu_0 p_0^2 \omega^4}{16 \pi c} \int_{S(r)} \left[ \sin^3 (\theta) d\theta \right] \end{aligned}

To calculate the integral of $\sin^3\theta$ , we can separate it into $\sin^2\theta * \sin\theta$ which equals $(1-\cos^2\theta)\sin\theta$ . Now we define $u = \cos\theta$ , so $du = - \sin(\theta) d\theta$ . Thus the integral of $\sin^3\theta$ is:

\begin{aligned} \int_{S(r)} (u^2 - 1) du & = \left[ \frac{1}{3}u^3 - u \right]_{S(r)} \\ & = \left[ \frac{1}{3}\cos^3\theta - \cos\theta \right]_0^\pi \\ & = \frac{2}{3} + \frac{2}{3} = \frac{4}{3} \end{aligned}

We finally have:

\begin{aligned} P & = \frac{\mu_0 p_0^2 \omega^4}{12 \pi c} \end{aligned}

Conclusion

Power vs Wavelength

We have $P \propto \frac{1}{\lambda^4}$ . We can see on the graph that in the visual spectrum $[400\:nm, 800\:nm]$ , the little wavelength (blues) get far more scattered than wavelength at the other end of the spectrum (reds), thus giving the sky its color.