Mathematics supporting the classical explanation of why the phase speed of light slows down in a medium
Edit : I tried to complete the Feynman calculation. Thank you for being kind to English, which is not my native language ! And sorry for possible miscalculation !
The problem is part of the Oseen extinction theorem. You'll find usueful reference in https://en.wikipedia.org/wiki/Ewald%E2%80%93Oseen_extinction_theorem In particular, the irremplace Born and Wolf. But it's complicated !
This is a simpler calculation. It is a macroscopic one. It is not perfect but I think it is interesting.
The field in a point is the sum of all the fields radiated by all the slices. But these radiated fields depend of the fied himself and so we are led to an integral equation.
In order not to be too long, I consider as known the field radiated by a plane of uniform polarization which oscillates sinusoidally.
A thin plate placed at $z=0$ and of length $dz$ is polarized : a volume $\text{d}\tau $ is equivalent to the dipole moment $\overrightarrow{\text{d}p}=\overrightarrow{P}\text{d}\tau =\overrightarrow{P}\text{d}z\text{d}S$. The polarization is supposed to vary sinusoidally in time and directed according to Ox: $\overrightarrow{P}=\overrightarrow{{{P}_{0}}}{{\operatorname{e}}^{j\omega t}}$with $\overrightarrow{{{P}_{0}}}={{P}_{0}}\overrightarrow{{{e}_{x}}}$ We work in complex representation but the complex amplitudes are not underlined to lighten the notations. All the vectors are along Ox.
The field radiated by this plate is ($k=\omega /c$):
$\text{d}{{E}_{x}}=-\frac{1}{2{{\varepsilon }_{0}}}(jk){{P}_{0}}dz{{\operatorname{e}}^{j\left( \omega t-kz \right)}}$ if $z>0$ and $\text{d}{{E}_{x}}=-\frac{1}{2{{\varepsilon }_{0}}}(jk){{P}_{0}}dz{{\operatorname{e}}^{j\left( \omega t+kz \right)}}$ if $z<0$
On can find the result in the Feynman's lectures : http://www.feynmanlectures.caltech.edu/I_30.html
I have simply adapted this result to a polarized medium.
Now consider a medium such that the polarization ${{P}_{x}}$ is related to the complex amplitude of the electric field by the relation ${{P}_{0}}={{\varepsilon }_{0}}\chi {{E}_{x}}$with the a priori complex quantity called susceptibility and dependent (a priori) of the frequency. The complex indices is linked to the susceptibility by ${{n}^{2}}=1+\chi $ . (This definition is easily adapted to a conductor with conductivity $\gamma $ )
The medium extends for z varying from 0 to infinity. For, $z<0$ we have vacuum. The medium is "lit" by an incident plane wave polarised along Ox, ${{E}_{i}}={{E}_{0}}{{\operatorname{e}}^{j}}^{(\omega t-kz)}$ with $k=\omega /c$. This wave will polarize the medium as in the previous question.
Each slice of the medium will then radiate. The total field at a point z, is the sum of the incident field and the fields radiated by the various slices.
By changing the origin for each slice, we have :
${{E}_{z'>z}}(z)=-\frac{1}{2{{\varepsilon }_{0}}}jk\int\limits_{z}^{+\infty }{{{P}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t+k(z-z') \right)}}}$
${{E}_{0<z'<z}}(z)=-\frac{1}{2{{\varepsilon }_{0}}}jk\int\limits_{0}^{z }{{{P}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t-k(z-z') \right)}}}$
Since ${{P}_{x}}(x')={{\varepsilon }_{0}}\chi {{E}_{x}}(x')$ :
${{E}_{x\,}}_{z'>z}(z)=-\frac{1}{2}jk\chi \int\limits_{z}^{+\infty }{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t+k(z-z') \right)}}}$
${{E}_{x\,0<z'<z}}(z)=-\frac{1}{2}jk\int\limits_{0}^{z }{\chi {{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t-k(z-z') \right)}}}$
At $z>0$ , the amplitude of the total field satisfies the integral equation: ${{E}_{x}}(z,t)={{E}_{0}}{{\operatorname{e}}^{j(\omega t-kz)}}\underbrace{-\frac{1}{2}jk\chi \int\limits_{z}^{+\infty }{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t+k(z-z') \right)}}}}_{{{E}_{x\,}}_{z'>z}(z)}\underbrace{-\frac{1}{2}jk\chi \int\limits_{0}^{z}{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{j\left( \omega t-k(z-z') \right)}}}}_{{{E}_{x\,}}_{z'<z}(z)}$
Since the medium is linear, we have ${{E}_{x}}(z,t)={{E}_{x}}(z){{\operatorname{e}}^{j(\omega t)}}$ and the field obey the integral equation, valid if $z>0$ :
${{E}_{x}}(z)={{E}_{0}}{{\operatorname{e}}^{-jkz}}-\frac{1}{2}jk\chi {{\operatorname{e}}^{jkz}}\int\limits_{z}^{+\infty }{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{-jkz'}}}-\frac{1}{2}jk\chi {{\operatorname{e}}^{-jkz}}\int\limits_{0}^{z}{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{+jkz'}}}$
To solve this integral equation, we can try a solution of the form ${{E}_{x}}(z)=C{{\operatorname{e}}^{\left( -j\beta z \right)}}$with C and $\beta $ constants a priori complex, and with $\operatorname{Im}(\beta )>0$ to avoid divergence. (We could also differentiate the equation twice).
$\int\limits_{0}^{z}{\text{C}{{\operatorname{e}}^{-j\beta z'}}\text{d}z'{{\operatorname{e}}^{+jkz'}}}=C\frac{1}{j(k-\beta )}({{\operatorname{e}}^{+j(k-\beta )z}}-1)$ $\int\limits_{z}^{\infty }{\text{C}{{\operatorname{e}}^{-j\beta z'}}\text{d}z'{{\operatorname{e}}^{-jkz'}}}=-C\frac{1}{j(k+\beta )}({{\operatorname{e}}^{-j(k+\beta )\infty }}-{{\operatorname{e}}^{-j(k+\beta )z}})=C\frac{1}{j(k+\beta )}{{\operatorname{e}}^{-j(k+\beta )z}}$
If we replace in the original equation:
$C{{\operatorname{e}}^{-j\beta z}}={{E}_{0}}{{\operatorname{e}}^{-jkz}}-\frac{1}{2}jk\chi {{\operatorname{e}}^{jkz}}C\frac{1}{j(k+\beta )}{{\operatorname{e}}^{-j(k+\beta )z}}-\frac{1}{2}jk\chi {{\operatorname{e}}^{-jkz}}C\frac{1}{j(k-\beta )}({{\operatorname{e}}^{+j(k-\beta )z}}-1)$
$C{{\operatorname{e}}^{-j\beta z}}={{E}_{0}}{{\operatorname{e}}^{-jkz}}+\frac{1}{2}jk\chi {{\operatorname{e}}^{-jkz}}\frac{C}{j(k-\beta )}-\frac{1}{2}jk\chi \frac{C}{j(k+\beta )}{{\operatorname{e}}^{-j\beta z}}-\frac{1}{2}jk\chi \frac{C}{j(k-\beta )}{{\operatorname{e}}^{-j\beta z}}$
We identify the term in ${{\operatorname{e}}^{-j\beta z}}$
$1=-\frac{1}{2}jk\chi \frac{1}{j(k+\beta )}-\frac{1}{2}jk\chi \frac{1}{j(k-\beta )}=-\chi \frac{{{k}^{2}}}{{{k}^{2}}-{{\beta }^{2}}}\to {{\beta }^{2}}={{n}^{2}}{{k}^{2}}={{n}^{2}}\frac{{{\omega }^{2}}}{{{c}^{2}}}$
So, we have $\beta =n\frac{\omega }{c}$
Identifying with the term in ${{\operatorname{e}}^{-jkz}}$ we find a second relation
$0={{E}_{0}}+\frac{1}{2}k\chi \frac{C}{(k-\beta )}$ . It gives $C=\frac{2(n-1)}{{{n}^{2}}-1}{{E}_{0}}=\frac{2}{n+1}{{E}_{0}}$
The field in the region $z>0$ is ${{E}_{t}}=\underbrace{\frac{2}{n+1}}_{t}{{E}_{0}}{{\operatorname{e}}^{j}}^{(\omega t-nkz)}$ and so we have the coefficient of transmission and the change of phase velocity by adding radiated fields with variable amplitudes, all propagating at velocity c.
But we can also compute the reflected waves in the $z<0$ region by adding all the radiated waves in this direction :
${{E}_{r}}(z)=-\frac{1}{2}jk\chi {{\operatorname{e}}^{jkz}}\int\limits_{0}^{+\infty }{{{E}_{x}}(z')\text{d}z'{{\operatorname{e}}^{-jkz'}}}=-{{\operatorname{e}}^{jkz}}\frac{{{n}^{2}}-1}{(n+1)}\frac{1}{n+1}{{E}_{0}}=-{{\operatorname{e}}^{jkz}}\frac{n-1}{n+1}{{E}_{0}}$
The reflected wave is ${{E}_{t}}(z,t)=\underbrace{\frac{1-n}{1+n}}_{r}{{E}_{0}}{{\operatorname{e}}^{j(\omega t+kz)}}$ and we have the coefficient of reflection without using any condition at the interface.
We can follow Feynman, chapter 31 of volume 1. The summary is :
A sinusoidal progressive plane wave in the vacuum ${{E}_{0}}{{e}^{i\omega (t-z/c)}}$ passes through an infinite plane of small thickness $\delta z$ and index $n$.
In the usual formalism of refraction index, it is delayed by $\delta t=(n-1)\delta z/c$ since it advances at the speed $c/n$ in the plate instead of $c$: time of travel in the plate $\delta z(n/c)$ instead of $\delta z(1/c)$.
The wave after the plate is $E={{E}_{0}}{{e}^{i(\omega (t-\delta t)-\omega z/c)}}={{e}^{i\omega (\delta t)}}{{E}_{0}}{{e}^{i(\omega (t-\delta t)-\omega z/c)}}=\underbrace{{{e}^{i\omega ((n-1)\delta z/c)}}}_{1-i\omega (n-1)\delta z/c}{{E}_{0}}{{e}^{i(\omega t-\omega z/c)}}$
Finally $E'={{E}_{0}}{{e}^{i(\omega t-z/c)}}\underbrace{-i\omega (n-1)\delta z/c{{E}_{0}}{{e}^{i(\omega t-z/c)}}}_{\text{Field radiated by the plane}}$
You just have to reverse the procedure: the sum of the incident wave and the radiated wave is indeed a delayed wave.
I recommend reading Feynman! (Sorry for my english)
The main idea really is that the changed phase velocity is a collective phenomenon that only manifests in macroscopic electric field, but is due to mutual microscopic interactions of the medium elements, while those interactions take place at unchanged vacuum light speed $c$.
The phase shift happens in the sense the farther the medium element is from the vacuum-medium interface, the larger the phase difference between the primary wave and the secondary wave generated at that element. The secondary waves cannot be generated in phase with the primary wave, because then they would, after some short path into the medium, screen the primary field completely and no macroscopic wave would penetrate deeper - we would be looking at total reflection of the wave, something that happens rather in metals, not in glass (at least not if wave direction is normal to the medium interface.
But the phase shift isn't something that immediately implies changed phase velocity, and it is not easy to find a detailed account of how that works.
Luckily, there is another and more clear argument based on macroscopic EM theory. In vacuum, the Maxwell equations imply the following wave equation for electric field:
$$ \frac{1}{c^2}\frac{\partial^2 E_x}{\partial t^2} - \frac{\partial^2E_x}{\partial z^2} =0. $$ If we restrict our attention to solutions which represent harmonic plane waves, all obey the universal relation: $$ \Omega / k = c $$ where $c$ is vacuum speed of light. I mention this because general definition of phase speed is $\Omega/k$. So, all such waves have phase speed $c$.
Similar wave equation can be derived for electric field in medium, but there is a difference - a new term. However, to derive this equation, we need to make some simplifying assumptions. These are:
1) second derivative of polarization of an element is proportional to instantaneous value of electric field there;
2) net macroscopic electric field inside medium can be expressed as harmonic plane wave, similarly to the vacuum case:
$$ \mathbf E = E_0 \mathbf e_x \sin(\Omega t - kz). $$
In such case, Maxwell's equations imply the following equation for electric field as a function of position $z$ and time $t$:
$$ \frac{1}{c^2}\frac{\partial^2 E_x}{\partial t^2} - \frac{\partial^2E_x}{\partial z^2} =-C\frac{\partial^2 E_x}{\partial t^2} $$ where $C$ is some constant quantifying how much the medium polarized in given macroscopic electric field.
Note the new term on the right-hand side, which was not present in the vacuum case. This terms accounts for the effect of polarized medium on the total electric field - the fact that polarized medium can generate its own, secondary radiation.
Solutions for $k$ as a function of $\Omega$ can be found by inserting the assumed form of electric field and some manipulation. It turns out that in general, due to the new term the ratio $\frac{\Omega}{k}$ does not have the vacuum value of $c$, but can be lower or higher, and is a function of $C$. This is described in terms of quantity $n$:
$$ \frac{\Omega}{k} = \frac{c}{n}. $$
In glass, usually $n>1$ so the waves are more densely packed along the path of propagation.
So, we can see how the changed phase speed is related to polarized medium and known laws. However, at the same time, we had to assume that waves of this simple kind are actually possible. This is not true in all cases, such as nonlinear media.
The key point is that the sum of two waves with different amplitude and phase but the same frequency gives you a wave of that frequency. This is shown in the picture where the black curve is $\sin(x)$, the red curve is $-.2 \sin(x+5)$ and the blue curve is the sum of the two. Same frequency, still a sine wave, but different phase.
Suppose there is an electric field $E=Ae^{i( \omega t-kx)}$. It falls on a molecule at (say) $x=0$. This is polarised by the field and oscillates and produces a field $Ere^{i \delta}$, where $r$ (which is small) depends on the polarisability of the molecule and the phase shift $\delta$ varies according to the difference between the applied $\omega$ and the natural resonant frequency of the molecule. So the combined field is
$Ae^{i \omega t }(1+re^{i \delta})$
The bracket can be written $(1 + r \cos\delta + i r \sin \delta)=Re^{i \phi}$
where $\tan \phi={r \sin \delta \over 1 + r \cos \delta}$ and $R^2=1+r^2+2 r \cos \delta$
So the effect of the polarisable molecule is just to shift the phase by a small $\phi$. (There is also a change in the amplitude but that's not at issue.)
If you have a series of such molecules then each has the same effect. The effect on the wave is multiplicative so the phases add. If there is a distance $a$ between molecules then in travelling a distance $x$ there are $x/a$ molecules giving a phase change due to polarisation of $\phi x/a$. The wave becomes $E=Ae^{i(\omega t - k x +\phi x/a)}=Ae^{i(\omega t -k' x)}$ with $k'=k-\phi/a$, i.e. the wavenumber (and thus the wavelength) change. The frequency is the same, so change in wavelength implies change in velocity.
This is the 'little bit of math' Orzel refers to.