diff --git a/src/4_sommerfeld_model.md b/src/4_sommerfeld_model.md
index 460700fbb9fd841e51254bdbab8e7e9cbb151d31..8fda87f5c6fc857e244e937304c642aed6cecb6c 100644
--- a/src/4_sommerfeld_model.md
+++ b/src/4_sommerfeld_model.md
@@ -29,23 +29,37 @@ _(based on chapter 4 of the book)_
     - calculate the electron contribution to the specific heat of a solid.
     - describe central terms such as the Fermi energy, Fermi temperature, and Fermi wavevector.
 
+## Electrons vs phonons
+
 > Two electrons are sitting on a bench. Another one approaches and asks: "May I join you guys?"
 > The first two immediately reply: "Who do you think we are? Bosons?"
 
-## Quantization of k-space and density of states in the free electron model
+Having learned the statistical properties of phonons and the [Debye model](2_debye_model.md), let us use these as a starting point for comparing with the electrons.
+Here is a table comparing the most important properties of electrons and phonons:
 
-Atoms in a metal provide conduction electrons from their outer shells (often s-shells). These electrons can be described as waves, analogous to phonons. The Hamiltonian of a free electron is:
+|   | Phonons | Electrons |
+| - | - | - |
+| Governed by | Wave equation | Schrödinger equation |
+|  | $ d² δ\mathbf{r} / dt² = -v²∇²δ\mathbf{r}$ | $ i ħdψ/dt = -ħ²∇²ψ/2m$ |
+| Dispersion relation | $ω=v \|\mathbf{k}\|$ | $ε = ħ²k²/2m$ |
+| Statistics | Bose-Einstein | Fermi-Dirac |
+| $n =$ | $1/[\exp(βE) - 1]$ | $1/[\exp(β[E - μ]) + 1]$ |
+| Number per $\mathbf{k}$ | 3 (polarization) | 2 (spin) |
+| Total number | temperature-dependent | as many as there are |
 
-$$
-\mathcal{H}=\frac{ {\bf p}^2}{2m}=-\frac{\hbar^2}{2m}\left( \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2} \right)\ \Rightarrow\ \varepsilon=\frac{\hbar^2}{2m}\left( k_x^2+k_y^2+k_z^2 \right)
-$$
+On the one hand, there are important differences, but on the other hand there are a lot of similarities.
 
-If we impose periodic boundary conditions, $\psi(x,y,z)=\psi(x+L,y+L,z+L)$, it follows that the allowed wavenumbers are given by
+* The solutions of the equation of motion are still plane waves $ψ ∝ \exp(i\mathbf{k}\mathbf{r})$.
+* The periodic boundary conditions work exactly the same: $k_{x,y,z}=…, \frac{-4\pi}{L}, \frac{-2\pi}{L}, 0, \frac{2\pi}{L}, \frac{4\pi}{L}, …$
+* The density of $k$-states per unit volume in the reciprocal space is the same: $(2π/L)³$.
+* The expression for the total energy of the system is *very similar*:  
+  $E_\textrm{total} = 2_s (2π/L)³∫ ε(\mathbf{k}) n_F[ε(\mathbf{k})]d³\mathbf{k}.$
 
-$$
-k_x=\frac{2\pi p}{L},\ k_y=\frac{2\pi q}{L},\ k_z=\frac{2\pi r}{L}\ \Rightarrow\ \varepsilon=\frac{2\pi^2\hbar^2}{mL^2}\left( p^2+q^2+r^2 \right).
-$$
-(this quantization of k-space is exactly analogous to that for phonons, see [lecture 2](2_debye_model.md) - Debye model)
+With this comparison we have everything to analyze the electron behavior.
+
+## Fermi sea
+
+For a warm-up let us imagine what happens in 1D, when we only have a few electron states available.
 
 ```python
 kf = 3;
@@ -67,7 +81,7 @@ ax.set_xlim(-3.75, 3.75);
 ax.set_ylim(0.0, 11);
 
 ax.set_xlabel(r"$k \enspace \left[ \frac{2 \pi}{L} \right]$");
-ax.set_ylabel(r"$\varepsilon$");
+ax.set_ylabel(r"$ε$");
 
 ax.set_xticklabels([""] + ks.tolist() + [""]);
 ax.set_yticks([]);
@@ -75,28 +89,31 @@ ax.set_yticks([]);
 draw_classic_axes(ax, xlabeloffset=.6);
 ```
 
-The difference with phonons is that electrons are _fermions_, implying that there are only 2 electron states (due to spin) per $k$-value.
+At $T=0$ the states with zero energy get completely filled up because the electrons cannot disappear.
 
 ![](figures/fermi_circle_periodic.svg)
 
-Our first goal is to compute the density of states. We start by expressing the total number of states $N$ as an integral over k-space. Assuming three dimensions and spherical symmetry (the dispersion in the free electron model is isotropic), we find
+## Density of states
+
+Our first goal is to compute the density of states, which we will do now using a more explicit algorithm.
+By definition density of states is the number of states per energy interval.
+Therefore if we compute the *total* number of states $N(ε)$ with energy lower than $ε$, then $g(ε) = dN(ε)/dε$.
 
+Let us apply this to the 3D case for a start.
+Assuming three dimensions and spherical symmetry (the dispersion in the free electron model is isotropic), we find
 $$
-N=2\left(\frac{L}{2\pi}\right)^3\int{\rm d}{\bf k}=2 \left(\frac{L}{2\pi}\right)^34\pi\int k^2{\rm d}k=\frac{V}{\pi^2}\int k^2{\rm d}k ,
+N=2\left(\frac{L}{2\pi}\right)^3\int{\rm d}{\bf k}=2 \left(\frac{L}{2\pi}\right)^34\pi\int k^2{\rm d}k=\frac{V}{\pi^2}\int k^2{\rm d}k,
 $$
+where the factor 2 is due to spin, and $\left(\frac{L}{2\pi}\right)^3$ is once again the density of points in k-space.
 
-where the factor 2 represents spin degeneracy, and $\left(\frac{L}{2\pi}\right)^3$ is the density of points in k-space.
-
-Using $k=\frac{\sqrt{2m\varepsilon}}{\hbar}$ and ${\rm d}k=\frac{1}{\hbar}\sqrt{\frac{m}{2\varepsilon}}{\rm d}\varepsilon$ we can rewrite this as:
-
+Using $k=\frac{\sqrt{2mε}}{\hbar}$ and ${\rm d}k=\frac{1}{\hbar}\sqrt{\frac{m}{2ε}}{\rm d}ε$ we can rewrite this as:
 $$
-N=\frac{V}{\pi^2}\int\frac{2m\varepsilon}{\hbar^3}\sqrt{\frac{m}{2\varepsilon}}{\rm d}\varepsilon=\frac{Vm^{3/2}}{\pi^2\hbar^3}\int\sqrt{2\varepsilon}\ {\rm d}\varepsilon
+N=\frac{V}{\pi^2}\int\frac{2mε}{\hbar^3}\sqrt{\frac{m}{2ε}}{\rm d}ε=\frac{Vm^{3/2}}{\pi^2\hbar^3}\int\sqrt{2ε}\ {\rm d}ε
 $$
 
 So we find for the density of states:
-
 $$
-g(\varepsilon)=\frac{ {\rm d}N}{ {\rm d}\varepsilon}=\frac{Vm^{3/2}\sqrt{2\varepsilon}}{\pi^2\hbar^3}\propto\sqrt{\varepsilon}
+g(ε)=\frac{ {\rm d}N}{ {\rm d}ε}=\frac{Vm^{3/2}\sqrt{2ε}}{\pi^2\hbar^3}\propto\sqrt{ε}
 $$
 
 ```python
@@ -105,38 +122,38 @@ fig, ax = pyplot.subplots()
 
 ax.plot(E, np.sqrt(E))
 
-ax.set_ylabel(r"$g(\varepsilon)$")
-ax.set_xlabel(r"$\varepsilon$")
+ax.set_ylabel(r"$g(ε)$")
+ax.set_xlabel(r"$ε$")
 draw_classic_axes(ax, xlabeloffset=.2)
 ```
 
 Similarly,
 
-- For 1D: $g(\varepsilon) = \frac{2 L}{\pi} \frac{ {\rm d}k}{ {\rm d}\varepsilon} \propto 1/\sqrt{\varepsilon}$
-- For 2D: $g(\varepsilon) = \frac{k L^2}{\pi} \frac{ {\rm d}k}{ {\rm d}\varepsilon} \propto \text{constant}$
+- For 1D: $g(ε) = \frac{2 L}{\pi} \frac{ {\rm d}k}{ {\rm d}ε} \propto 1/\sqrt{ε}$
+- For 2D: $g(ε) = \frac{k L^2}{\pi} \frac{ {\rm d}k}{ {\rm d}ε} \propto \text{constant}$
 
 
 Given the number of electrons in a system, we can now fill up these states starting from the lowest energy until we run out of electrons, at which point we reach the _Fermi energy_.
 
 ## Fermi energy, Fermi wavevector, Fermi wavelength
 
-At $T=0$, the total number of electrons is given by the integral over the density of states up to the _Fermi energy_ $\varepsilon_{\rm F}$:
+At $T=0$, the total number of electrons is given by the integral over the density of states up to the _Fermi energy_ $ε_{\rm F}$:
 
 $$
-N=\int_0^{\varepsilon_{\rm F}}g(\varepsilon){\rm d}\varepsilon \overset{\mathrm{3D}}{=} \frac{V}{3\pi^2\hbar^3}(2m\varepsilon_{\rm F})^{3/2}.
+N=\int_0^{ε_{\rm F}}g(ε){\rm d}ε \overset{\mathrm{3D}}{=} \frac{V}{3\pi^2\hbar^3}(2mε_{\rm F})^{3/2}.
 $$
 Note that $N$ now denotes the total number of electrons in the system.
 
-Alternatively, we can express $N$ as an integral over k-space up to the _Fermi wavenumber_, which is the wavenumber associated with the Fermi energy $k_{\rm F}=\sqrt{2m\varepsilon_{\rm F}}/\hbar$
+Alternatively, we can express $N$ as an integral over k-space up to the _Fermi wavenumber_, which is the wavenumber associated with the Fermi energy $k_{\rm F}=\sqrt{2mε_{\rm F}}/\hbar$
 
 $$
 N \overset{\mathrm{3D}}{=} 2\frac{L^3}{(2\pi)^3}\int_0^{k_{\rm F}} 4\pi k^2{\rm d}k= 2\frac{V}{(2\pi)^3} \frac{4}{3}\pi k_{\rm F}^3
 $$
 
-These equations allow us to relate $\varepsilon_{\rm F}$ and $k_{\rm F}$ to the electron density $N/V$:
+These equations allow us to relate $ε_{\rm F}$ and $k_{\rm F}$ to the electron density $N/V$:
 
 $$
-\varepsilon_{\rm F}=\frac{\hbar^2}{2m}\left( 3\pi^2\frac{N}{V} \right)^{2/3}, {\rm and} \quad k_{\rm F}=\left( 3\pi^2\frac{N}{V} \right)^{1/3}.
+ε_{\rm F}=\frac{\hbar^2}{2m}\left( 3\pi^2\frac{N}{V} \right)^{2/3}, {\rm and} \quad k_{\rm F}=\left( 3\pi^2\frac{N}{V} \right)^{1/3}.
 $$
 
 From the last equation it follows that the _Fermi wavelength_ $\lambda_{\rm F}\equiv 2\pi/k_{\rm F}$ is on the order of the atomic spacing for typical free electron densities in metals.
@@ -164,10 +181,10 @@ ax.axhline(kf*kf, linestyle="dotted", color='k');
 ax.set_xticks([kf]);
 ax.set_yticks([kf*kf + 0.4]);
 ax.set_xticklabels([r"$k_F$"]);
-ax.set_yticklabels([r"$\varepsilon_F$"]);
+ax.set_yticklabels([r"$ε_F$"]);
 
 ax.set_xlabel(r"$k$");
-ax.set_ylabel(r"$\varepsilon$");
+ax.set_ylabel(r"$ε$");
 
 ax.set_xlim(-kf*extrapol, kf*extrapol);
 ax.set_ylim(0.0, kf*kf*extrapol);
@@ -176,7 +193,7 @@ draw_classic_axes(ax, xlabeloffset=.6);
 
 The bold line represents all filled states at $T=0$. This is called the _Fermi sea_.
 
-New concept: _Fermi surface_ = all points in k-space with $\varepsilon=\varepsilon_{\rm F}$. For free electrons in 3D, the Fermi surface is the surface of a sphere.
+New concept: _Fermi surface_ = all points in k-space with $ε=ε_{\rm F}$. For free electrons in 3D, the Fermi surface is the surface of a sphere.
 
 
 ![](figures/transport.svg)
@@ -185,16 +202,16 @@ The orange circle represents the Fermi surface at finite current $\rightarrow$ t
 
 
 ## Finite temperature, heat capacity
-We now extend our discussion to $T>0$ by including a temperature dependent occupation function $n_F(\varepsilon,T)$ into our expression for the total number of electrons:
+We now extend our discussion to $T>0$ by including a temperature dependent occupation function $n_F(ε,T)$ into our expression for the total number of electrons:
 
 $$
-N=\int_0^\infty n_F(\varepsilon,T)g(\varepsilon){\rm d}\varepsilon,
+N=\int_0^\infty n_F(ε,T)g(ε){\rm d}ε,
 $$
 
 where the probability for a certain electron state to be occupied is given by the Fermi-Dirac distribution
 
 $$
-n_F(\varepsilon,T)=\frac{1}{ {\rm e}^{(\varepsilon-\mu)/k_{\rm B}T}+1}
+n_F(ε,T)=\frac{1}{ {\rm e}^{(ε-\mu)/k_{\rm B}T}+1}
 $$
 
 ```python
@@ -206,8 +223,8 @@ beta = 20
 ax.plot(xvals, xvals < mu, ls='dashed', label='$T=0$')
 ax.plot(xvals, 1/(np.exp(beta * (xvals-mu)) + 1),
         ls='solid', label='$T>0$')
-ax.set_xlabel(r'$\varepsilon$')
-ax.set_ylabel(r'$n_F(\varepsilon, T)$')
+ax.set_xlabel(r'$ε$')
+ax.set_ylabel(r'$n_F(ε, T)$')
 ax.set_yticks([0, 1])
 ax.set_yticklabels(['$0$', '$1$'])
 ax.set_xticks([mu])
@@ -217,7 +234,7 @@ ax.legend()
 draw_classic_axes(ax)
 pyplot.tight_layout()
 ```
-where the chemical potential $\mu=\varepsilon_{\rm F}$ if $T=0$. Typically $\varepsilon_{\rm F}/k_{\rm B}$~70 000 K (~7 eV), whereas room temperature is only 300 K (~30 meV). Therefore, thermal smearing occurs only very close to the Fermi energy.
+where the chemical potential $\mu=ε_{\rm F}$ if $T=0$. Typically $ε_{\rm F}/k_{\rm B}$~70 000 K (~7 eV), whereas room temperature is only 300 K (~30 meV). Therefore, thermal smearing occurs only very close to the Fermi energy.
 
 Having included the temperature dependence, we can now calculate the electronic contribution to the heat capacity $C_{\rm V,e}$.
 
@@ -225,7 +242,7 @@ Having included the temperature dependence, we can now calculate the electronic
 E = np.linspace(0, 2, 500)
 fig, ax = pyplot.subplots()
 ax.plot(E, np.sqrt(E), linestyle='dashed')
-ax.text(1.7, 1.4, r'$g(\varepsilon)\propto \sqrt{\varepsilon}$', ha='center')
+ax.text(1.7, 1.4, r'$g(ε)\propto \sqrt{ε}$', ha='center')
 ax.fill_between(E, np.sqrt(E) * (E < 1), alpha=.3)
 
 n = np.sqrt(E) / (1 + np.exp(20*(E-1)))
@@ -238,29 +255,29 @@ ax.text(1-w/2, 1.1, r'$\sim k_BT$', ha='center')
 ax.plot([1-w, 1+w], [1, 0], c='k', linestyle='dashed')
 ax.annotate(s='', xy=(1, 0), xytext=(1, 1),
             arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-ax.text(1.2, .7, r'$g(\varepsilon_F)$', ha='center')
+ax.text(1.2, .7, r'$g(ε_F)$', ha='center')
 ax.set_xticks([1])
-ax.set_xticklabels([r'$\varepsilon_F$'])
+ax.set_xticklabels([r'$ε_F$'])
 
-ax.set_ylabel(r"$g(\varepsilon)$")
-ax.set_xlabel(r"$\varepsilon$")
+ax.set_ylabel(r"$g(ε)$")
+ax.set_xlabel(r"$ε$")
 draw_classic_axes(ax, xlabeloffset=.2)
 ```
 
-We will estimate $C_{\rm V,e}$ and derive its scaling with temperature using the triangle method depicted in the figure. A finite temperature causes electrons in the top triangle to be excited to the bottom triangle. Because the base of this triangle scales with $k_{\rm B}T$ and its height with $ g(\varepsilon_{\rm F})$, it follows that the number of excited electrons $N_{\rm exc} \approx g(\varepsilon_{\rm F})k_{\rm B}T$ (neglecting pre-factors of order 1).
+We will estimate $C_{\rm V,e}$ and derive its scaling with temperature using the triangle method depicted in the figure. A finite temperature causes electrons in the top triangle to be excited to the bottom triangle. Because the base of this triangle scales with $k_{\rm B}T$ and its height with $ g(ε_{\rm F})$, it follows that the number of excited electrons $N_{\rm exc} \approx g(ε_{\rm F})k_{\rm B}T$ (neglecting pre-factors of order 1).
 
 These electrons gain $k_{\rm B}T$ of energy, so the total extra energy is
 
 $$
-E(T)-E(0)=N_{\rm exc}k_{\rm B}T\approx g(\varepsilon_{\rm F})k_{\rm B}^2T^2.
+E(T)-E(0)=N_{\rm exc}k_{\rm B}T\approx g(ε_{\rm F})k_{\rm B}^2T^2.
 $$
 
 Therefore, the heat capacity is given by
 
 $$
-C_{\rm V,e}=\frac{ {\rm d}E}{ {\rm d}T} \approx 2 g(\varepsilon_{\rm F})k_{\rm B}^2T=\ ...\ =3 Nk_{\rm B}\frac{T}{T_{\rm F}}\propto T,
+C_{\rm V,e}=\frac{ {\rm d}E}{ {\rm d}T} \approx 2 g(ε_{\rm F})k_{\rm B}^2T=\ ...\ =3 Nk_{\rm B}\frac{T}{T_{\rm F}}\propto T,
 $$
-where we used $N=\frac{2}{3}\varepsilon_{\rm F}g(\varepsilon_{\rm F})$ and we defined the _Fermi temperature_ $T_{\rm F}=\frac{\varepsilon_{\rm F}}{k_{\rm B}}$.
+where we used $N=\frac{2}{3}ε_{\rm F}g(ε_{\rm F})$ and we defined the _Fermi temperature_ $T_{\rm F}=\frac{ε_{\rm F}}{k_{\rm B}}$.
 
 How does $C_{\rm V,e}$ relate to the phonon contribution $C_{\rm V,p}$?
 
@@ -275,19 +292,19 @@ Behavior of $C_{\rm V}$ can be very quickly memorized or understood using the fo
 
 #### Example 1: electrons
 
-$g(\varepsilon_{\rm F})$ roughly constant ⇒ total energy in the thermal state is $T \times [T\times g(\varepsilon_{\rm F})]$ ⇒ $C_{\rm V} \propto T$.
+$g(ε_{\rm F})$ roughly constant ⇒ total energy in the thermal state is $T \times [T\times g(ε_{\rm F})]$ ⇒ $C_{\rm V} \propto T$.
 
 #### Example 2: graphene with $E_F=0$ (midterm 2018)
 
-$g(\varepsilon) \propto \varepsilon$ ⇒ total energy is $T \times T^2$ ⇒ $C_{\rm V} \propto T^2$.
+$g(ε) \propto ε$ ⇒ total energy is $T \times T^2$ ⇒ $C_{\rm V} \propto T^2$.
 
 #### Example 3: phonons in 3D at low temperatures.
 
-$g(\varepsilon) \propto \varepsilon^2$ ⇒ total energy is $T \times T^3$ ⇒ $C_{\rm V} \propto T^3$.
+$g(ε) \propto ε^2$ ⇒ total energy is $T \times T^3$ ⇒ $C_{\rm V} \propto T^3$.
 
 
 ## Conclusions
-  1. The Sommerfeld free electron model treats electrons as waves with dispersion $\varepsilon=\frac{\hbar^2k^2}{2m}$.
+  1. The Sommerfeld free electron model treats electrons as waves with dispersion $ε=\frac{\hbar^2k^2}{2m}$.
   2. The density of states (DOS) can be derived from the dispersion relation. This procedure is general, and analogous to e.g. that for phonons (see lecture 2 - Debye model).
   3. The Fermi-Dirac distribution describes the probability an electron state is occupied.
   4. The free-electron heat capacity is linear with $T$.