use inline code execution

docs/lecture_1.md → src/lecture_1.md

+```python
+from matplotlib import pyplot
+
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy.integrate import quad
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+```
 
 # Lecture 1 – Phonons and specific Heat
 _Based on chapter 2 of the book_
@@ -25,7 +36,31 @@ $$
 n(\omega,T)=\frac{1}{ {\rm e}^{\hbar\omega/k_{\rm B}T}-1}\Rightarrow\bar{\varepsilon}=\frac{1}{2}\hbar\omega+\frac{\hbar\omega}{ {\rm e}^{\hbar\omega/k_{\rm B}T}-1}
 $$
 
-![](figures/bose_einstein.svg)
+```python
+
+fig, (ax, ax2) = pyplot.subplots(ncols=2, figsize=(10, 5))
+omega = np.linspace(0.1, 2)
+ax.plot(omega, 1/(np.exp(omega) - 1))
+ax.set_ylim(0, top=3)
+ax.set_xlim(left=0)
+ax.set_xlabel(r'$\hbar \omega$')
+ax.set_xticks([0, 1])
+ax.set_xticklabels(['$0$', '$k_B T$'])
+ax.set_ylabel('$n$')
+ax.set_yticks([1, 2])
+ax.set_yticklabels(['$1$', '$2$'])
+draw_classic_axes(ax, xlabeloffset=.2)
+T = np.linspace(0.01, 2)
+ax2.plot(T, 1/2 + 1/(np.exp(1/T)-1))
+ax2.set_ylim(bottom=0)
+ax2.set_xlabel('$k_B T$')
+ax2.set_xticks([0, 1])
+ax2.set_xticklabels(['$0$', r'$\hbar \omega$'])
+ax2.set_ylabel(r"$\bar\varepsilon$")
+ax2.set_yticks([1/2])
+ax2.set_yticklabels([r'$\hbar\omega/2$'])
+draw_classic_axes(ax2, xlabeloffset=.15)
+```
 
 The term $\frac{1}{2}\hbar\omega$ is the _zero point energy_, which follows from the uncertainty principle.
 
@@ -40,7 +75,31 @@ C = \frac{\partial\bar{\varepsilon}}{\partial T}
 \end{multline}
 $$
 
-![](figures/zeropoint.svg)
+```python
+# Data from Einstein's paper
+T = [222.4, 262.4, 283.7, 306.4, 331.3, 358.5, 413.0, 479.2, 520.0, 879.7, 1079.7, 1258.0],
+c = [0.384, 0.578, 0.683, 0.798, 0.928, 1.069, 1.343, 1.656, 1.833, 2.671, 2.720, 2.781]
+
+def c_einstein(T, T_E=1):
+    x = T_E / T
+    return 3 * x**2 * np.exp(x) / (np.exp(x) - 1)**2
+
+T = np.linspace(0.01, 1.5, 500)
+fig, ax = pyplot.subplots()
+
+ax.plot(T, c_einstein(T)/3)
+ax.fill_between(T, c_einstein(T)/3, 1, alpha=0.5)
+
+ax.set_ylim(bottom=0, top=1.2)
+ax.set_xlabel('$T$')
+ax.set_ylabel(r'$\omega$')
+ax.set_xticks([1])
+ax.set_xticklabels([r'$\hbar \omega/k_B$'])
+ax.set_yticks([1])
+ax.set_yticklabels(['$k_B$'])
+pyplot.hlines([1], 0, 1.5, linestyles='dashed')
+draw_classic_axes(ax)
+```
 
 The dashed line signifies the classical value, $k_{\rm B}$. Shaded area $=\frac{1}{2}\hbar\omega$, the zero point energy that cannot be removed through cooling.
 
@@ -126,4 +185,29 @@ $$
 
 where $x=\frac{\hbar\omega}{k_{\rm B}T}$ and $\Theta_{\rm D}\equiv\frac{\hbar\omega_{\rm D}}{k_{\rm B}}$, the _Debye temperature_.
 
-![](figures/debye.svg)
+```python
+def integrand(y):
+    return y**4 * np.exp(y) / (np.exp(y) - 1)**2
+
+@np.vectorize
+def c_debye(T, T_D=1):
+    x = T / T_D
+    return 9 * x**3 * quad(integrand, 0, 1/x)[0]
+
+T = np.linspace(0.01, 1.5, 500)
+fig, ax = pyplot.subplots()
+
+ax.plot(T, c_einstein(T), label="Einstein model")
+ax.plot(T, c_debye(T), label="Debye model")
+
+ax.set_ylim(bottom=0, top=3.4)
+ax.set_xlabel('$T$')
+ax.set_ylabel(r'$\omega$')
+ax.set_xticks([1])
+ax.set_xticklabels([r'$\Theta_D$'])
+ax.set_yticks([3])
+ax.set_yticklabels(['$3Nk_B$'])
+ax.legend(loc='lower right')
+pyplot.hlines([3], 0, 1.5, linestyles='dashed')
+draw_classic_axes(ax, xlabeloffset=0.3)
+```

docs/lecture_2.md → src/lecture_2.md

+```python
+from matplotlib import pyplot
+
+import numpy as np
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+```
 
 # Lecture 2 – Free electron model
 _Based on chapters 3 and 4 of the book_
@@ -51,7 +60,7 @@ How fast do electrons travel through a copper wire? Let's take $E$ = 1 volt/m, $
 $\rightarrow v=\mu E=\frac{e\tau}{m_{\rm e}}E=\frac{10^{-19}\times 2.5\times 10^{-14}}{10^{-30}}=2.5\times10^{-3}=2.5$ mm/s ! (= 50 $\mu$m @ 50 Hz AC)
 
 ### Hall effect
-Consider a conductive wire in a magnetic field ${\bf B} \rightarrow$ electrons are deflected in a direction perpendicular to ${\bf B}$ and ${\bf j}$. 
+Consider a conductive wire in a magnetic field ${\bf B} \rightarrow$ electrons are deflected in a direction perpendicular to ${\bf B}$ and ${\bf j}$.
 
 ![](figures/hall_effect.svg)
 
@@ -94,7 +103,7 @@ Comparable to phonons, but: electrons are _fermions_.
 
 - Only 2 (due to spin) allowed per $k$-value
 - Fill up from the lowest energy until you run out of electrons
-   
+
 $\rightarrow$ Calculate when you are out of electrons $\rightarrow$ _Fermi energy_.
 
 ![](figures/fermi_circle_periodic.svg)
@@ -117,7 +126,16 @@ $$
 g(\varepsilon)=\frac{ {\rm d}N}{ {\rm d}\varepsilon}=\frac{Vm^{3/2}\sqrt{2\varepsilon}}{\pi^2\hbar^3}\propto\sqrt{\varepsilon}
 $$
 
-![](figures/sqrt.svg)
+```python
+E = np.linspace(0, 2, 500)
+fig, ax = pyplot.subplots()
+
+ax.plot(E, np.sqrt(E))
+
+ax.set_ylabel(r"$g(\varepsilon)$")
+ax.set_xlabel(r"$\varepsilon$")
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 Similarly,
 
@@ -133,7 +151,7 @@ $$
 with $\varepsilon_{\rm F}$ the _Fermi energy_ = highest filled energy at $T=0$.
 
 $$
-\varepsilon_{\rm F}=\frac{\hbar^2}{2m}\left( 3\pi^2\frac{N}{V} \right)^{2/3}\equiv \frac{\hbar^2 k_{\rm F}^2}{2m},\ 
+\varepsilon_{\rm F}=\frac{\hbar^2}{2m}\left( 3\pi^2\frac{N}{V} \right)^{2/3}\equiv \frac{\hbar^2 k_{\rm F}^2}{2m},\
 k_{\rm F}=\left( 3\pi^2\frac{N}{V} \right)^{1/3}
 $$
 
@@ -155,7 +173,27 @@ $$
 f(\varepsilon,T)=\frac{1}{ {\rm e}^{(\varepsilon-\mu)/k_{\rm B}T}+1}
 $$
 
-![](figures/fermi_dirac.svg)
+```python
+fig = pyplot.figure()
+ax = fig.add_subplot(1,1,1)
+xvals = np.linspace(0, 2, 200)
+mu = .75
+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'$f(\varepsilon, T)$')
+ax.set_yticks([0, 1])
+ax.set_yticklabels(['$0$', '$1$'])
+ax.set_xticks([mu])
+ax.set_xticklabels([r'$\mu$'])
+ax.set_ylim(-.1, 1.1)
+ax.legend()
+draw_classic_axes(ax)
+pyplot.tight_layout()
+```
+
 
 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) $\rightarrow$ thermal smearing occurs only very close to Fermi surface.
 
@@ -167,7 +205,31 @@ $$
 
 We can use this to calculate the electronic contribution to the heat capacity.
 
-![](figures/FD_DOS.svg)
+```python
+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.fill_between(E, np.sqrt(E) * (E < 1), alpha=.3)
+
+n = np.sqrt(E) / (1 + np.exp(20*(E-1)))
+ax.plot(E, n)
+ax.fill_between(E, n, alpha=.5)
+w = .17
+ax.annotate(s='', xy=(1, 1), xytext=(1-w, 1),
+            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+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.set_xticks([1])
+ax.set_xticklabels([r'$\varepsilon_F$'])
+
+ax.set_ylabel(r"$g(\varepsilon)$")
+ax.set_xlabel(r"$\varepsilon$")
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 Electrons in the top triangle are being excited to the bottom triangle due to temperature increase. Number of excited electrons $\approx\frac{1}{2}g(\varepsilon_{\rm F})k_{\rm B}T=n_{\rm exc}$. Total extra energy $E(T)-E(0)=n_{\rm exc}k_{\rm B}T=\frac{1}{2}g(\varepsilon_{\rm F})k_{\rm B}^2T^2$.
 

docs/lecture_3.md → src/lecture_3.md

+```python
+from matplotlib import pyplot
+
+import numpy as np
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+
+pi = np.pi
+```
 
 # Atoms and bonds
 
@@ -246,7 +257,7 @@ Same as a single covalent bond, only more atoms in a line. Considering 3 atoms:
 $$
 E \begin{pmatrix}
 c_1 \\ c_2 \\ c_3
-\end{pmatrix} = 
+\end{pmatrix} =
 \begin{pmatrix}
 E_0 & t & 0 \\
 t & E_0 & t \\
@@ -263,22 +274,45 @@ Diagonalizing large matrices is unwieldy, but let's try and check it numerically
 
 Eigenfrequencies of 3 atoms: `[0.0 1.0 1.732050]`
 
-![](figures/3_phonons.svg)
+```python
+def DOS_finite_phonon_chain(n):
+    rhs = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
+    rhs[0, 0] -= 1
+    rhs[-1, -1] -= 1
+    pyplot.figure()
+    pyplot.hist(np.sqrt(np.abs(np.linalg.eigvalsh(rhs))), bins=30)
+    pyplot.xlabel("$\omega$")
+    pyplot.ylabel("Number of levels")
+
+DOS_finite_phonon_chain(3)
+```
 
 Energies of 3 orbitals: `[-1.41421356 0.0  1.41421356]`
 
-![](figures/3_electrons.svg)
+```python
+def DOS_finite_electron_chain(n):
+    rhs = - np.eye(n, k=1) - np.eye(n, k=-1)
+    pyplot.figure()
+    pyplot.hist(np.linalg.eigvalsh(rhs), bins=30)
+    pyplot.xlabel("$E$")
+    pyplot.ylabel("Number of levels")
+
+DOS_finite_electron_chain(3)
+```
 
 ### From 3 atoms to 300
 
 Frequencies of many phonons:
 
-![](figures/300_phonons.svg)
+```python
+DOS_finite_phonon_chain(300)
+```
 
 Energies of many electrons:
 
-![](figures/300_electrons.svg)
-
+```python
+DOS_finite_electron_chain(300)
+```
 ## Summary
 
 * Electrons in atoms occupy "shells" in the energetically favorable order

docs/lecture_4.md → src/lecture_4.md

+```python
+import matplotlib
+from matplotlib import pyplot
+
+import numpy as np
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+
+pi = np.pi
+```
+
 # Electrons and phonons in 1D
 
 Last lecture:
@@ -10,7 +23,7 @@ This lecture:
 * Infinitely many atoms
 * Main idea: use *periodicity* in space, similar to periodicity in time
 
-## Equations of motion 
+## Equations of motion
 
 ### Phonons
 
@@ -53,7 +66,24 @@ In that case $c_n/c_0 = \exp(2 \pi n p a/L) = \exp(2 \pi n p/N)$, and changing $
 
 We can see that the solutions with different $k$ (or different $p$) are identical by plotting them:
 
-![](figures/phonons4.svg)
+```python
+x = np.linspace(-.2, 2.8, 500)
+fig, ax = pyplot.subplots()
+ax.plot(x, np.cos(pi*(x - .3)), label=r'$k=\pi/a$')
+ax.plot(x, np.cos(3*pi*(x - .3)), label=r'$k=3\pi/a$')
+ax.plot(x, np.cos(5*pi*(x - .3)), label=r'$k=5\pi/a$')
+sites = np.arange(3) + .3
+ax.scatter(sites, np.cos(pi*(sites - .3)), c='k', s=64, zorder=5)
+ax.set_xlabel('$x$')
+ax.set_ylabel('$u_n$')
+ax.set_xlim((-.1, 3.2))
+ax.set_ylim((-1.3, 1.3))
+ax.legend(loc='lower right')
+draw_classic_axes(ax)
+ax.annotate(s='', xy=(.3, -1.1), xytext=(1.3, -1.1),
+            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+ax.text(.3 + .5, -1.25, '$a$', ha='center');
+```
 
 Let's count how many different solutions we expect to find. We have a system with $N$ degrees of freedom (either $u_n$ or $c_n$), and therefore we need $N$ normal modes (or eigenstates).
 
@@ -73,7 +103,27 @@ $$\omega = \sqrt{\frac{2\kappa}{m}}\sqrt{1-\cos(ka)} = 2\sqrt{\frac{\kappa}{m}}|
 
 So we arrive to the dispersion relation
 
-![](figures/phonons3.svg)
+```python
+k = np.linspace(-2*pi, 6*pi, 500)
+fig, ax = pyplot.subplots()
+
+pyplot.plot(k, np.abs(np.sin(k/2)))
+
+ax.set_ylim(bottom=0, top=1.2)
+ax.set_xlabel('$ka$')
+ax.set_ylabel(r'$\omega$')
+pyplot.xticks(list(2*pi*np.arange(-1, 4)) + [-pi, pi],
+              [r'$-2\pi$', '$0$', r'$2\pi$', r'$4\pi$', r'$6\pi$',
+               r'$-\pi$', r'$\pi$'])
+pyplot.yticks([1], [r'$2\sqrt{\frac{\kappa}{m}}$'])
+pyplot.vlines([-pi, pi], 0, 1.1, linestyles='dashed')
+pyplot.hlines([1], .1*pi, 1.3*pi, linestyles='dashed')
+draw_classic_axes(ax)
+ax.annotate(s='', xy=(-pi, -.15), xytext=(pi, -.15),
+            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+ax.text(0, -.25, '1st Brillouin zone', ha='center')
+ax.set_ylim(bottom=-.7);
+```
 
 Here $k_p = 2\pi p / L$ for $0 ≤ p < N$ are *all the normal modes* of the crystal, and therefore we can rederive a better Debye model!
 
@@ -96,7 +146,14 @@ E = E_0 + 2t\cos(ka),
 $$
 so we come to the dispersion relation:
 
-![png](figures/tight_binding_1d.svg)
+```python
+pyplot.figure()
+k = np.linspace(-pi, pi, 300)
+pyplot.plot(k, -np.cos(k))
+pyplot.xlabel('$ka$'); pyplot.ylabel('$E$')
+pyplot.xticks([-pi, 0, pi], [r'$-\pi$', 0, r'$\pi$'])
+pyplot.yticks([-1, 0, 1], ['$E_0+2t$', '$E_0$', '$E_0-2t$']);
+```
 
 
 Usually electron dispersion has multiple options for $E(k)$, each called a *band*. The complete dispersion relation is also called a *band structure*.
@@ -147,8 +204,17 @@ Substituting, we get
 
 $$m_{eff} = \hbar^2\left(\frac{d^2 E(k)}{dk^2}\right)^{-1}$$
 
-![](figures/meff_1d_tb.svg)
-
+```python
+pyplot.figure(figsize=(8, 5))
+k = np.linspace(-pi, pi, 300)
+meff = 1/np.cos(k)
+color = list(matplotlib.rcParams['axes.prop_cycle'])[0]['color']
+pyplot.plot(k[meff > 0], meff[meff > 0], c=color)
+pyplot.plot(k[meff < 0], meff[meff < 0], c=color)
+pyplot.ylim(-5, 5)
+pyplot.xlabel('$ka$'); pyplot.ylabel('$m_{eff}$')
+pyplot.xticks([-pi, 0, pi], [r'$-\pi$', 0, r'$\pi$']);
+```
 
 ### density of states
 
@@ -251,7 +317,30 @@ $$\omega^2=\frac{\kappa(M+m)}{Mm}\pm \kappa\left\{\left(\frac{M+m}{Mm}\right)^2-
 
 Looking at the eigenvectors we see that for every $k$ there are now two values of $\omega$: one corresponding to in-phase motion (–) and anti-phase (+).
 
-![](figures/phonons6.svg)
+```python
+def dispersion_2m(k, kappa=1, M=1.4, m=1, acoustic=True):
+    Mm = M*m
+    m_harm = (M + m) / Mm
+    root = kappa * np.sqrt(m_harm**2 - 4*np.sin(k/2)**2 / Mm)
+    if acoustic:
+        root *= -1
+    return np.sqrt(kappa*m_harm + root)
+
+# TODO: Add panels with eigenvectors
+k = np.linspace(-2*pi, 6*pi, 300)
+fig, ax = pyplot.subplots()
+ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
+ax.plot(k, dispersion_2m(k), label='acoustic')
+ax.set_xlabel('$ka$')
+ax.set_ylabel(r'$\omega$')
+pyplot.xticks([-pi, 0, pi], [r'$-\pi$', '$0$', r'$\pi$'])
+pyplot.yticks([], [])
+pyplot.vlines([-pi, pi], 0, 2.2, linestyles='dashed')
+pyplot.legend()
+pyplot.xlim(-1.75*pi, 3.5*pi)
+pyplot.ylim(bottom=0)
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 'A' is called the _acoustic branch_, 'B' is called the _optical branch_.
 
@@ -261,9 +350,46 @@ The quantity $\frac{ {\rm d}k}{ {\rm d}\omega}$ can be calculated from the dispe
 
 When plotted, the DOS may look something like this:
 
-![](figures/phonons8.svg)
-
-Note that normally $g(\omega)$ is plotted along the vertical axis and $\omega$ along the horizontal axis – the plot above is just to demonstrate the relation between the dispersion and the DOS. The singularities in $g(\omega)$ at the bottom and top of each branch are called _van Hove singularities_. 
+```python
+matplotlib.rcParams['font.size'] = 24
+k = np.linspace(-.25*pi, 1.5*pi, 300)
+k_dos = np.linspace(0, pi, 20)
+fig, (ax, ax2) = pyplot.subplots(ncols=2, sharey=True, figsize=(10, 5))
+ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
+ax.plot(k, dispersion_2m(k), label='acoustic')
+ax.vlines(k_dos, 0, dispersion_2m(k_dos, acoustic=False),
+          colors=(0.5, 0.5, 0.5, 0.5))
+ax.hlines(
+    np.hstack((dispersion_2m(k_dos, acoustic=False), dispersion_2m(k_dos))),
+    np.hstack((k_dos, k_dos)),
+    1.8*pi,
+    colors=(0.5, 0.5, 0.5, 0.5)
+)
+ax.set_xlabel('$ka$')
+ax.set_ylabel(r'$\omega$')
+ax.set_xticks([0, pi])
+ax.set_xticklabels(['$0$', r'$\pi$'])
+ax.set_yticks([], [])
+ax.set_xlim(-pi/4, 2*pi)
+ax.set_ylim((0, dispersion_2m(0, acoustic=False) + .2))
+draw_classic_axes(ax, xlabeloffset=.2)
+
+k = np.linspace(0, pi, 1000)
+omegas = np.hstack((
+    dispersion_2m(k, acoustic=False), dispersion_2m(k)
+))
+ax2.hist(omegas, orientation='horizontal', bins=75)
+ax2.set_xlabel(r'$g(\omega)$')
+ax2.set_ylabel(r'$\omega$')
+
+# Truncate the singularity in the DOS
+max_x = ax2.get_xlim()[1]
+ax2.set_xlim((0, max_x/2))
+draw_classic_axes(ax2, xlabeloffset=.1)
+matplotlib.rcParams['font.size'] = 16
+```
+
+Note that normally $g(\omega)$ is plotted along the vertical axis and $\omega$ along the horizontal axis – the plot above is just to demonstrate the relation between the dispersion and the DOS. The singularities in $g(\omega)$ at the bottom and top of each branch are called _van Hove singularities_.
 
 ### Consistency check with 1 atom per cell
 
@@ -275,7 +401,25 @@ For $M\rightarrow m$, the dispersion diagram should come back to the diagram obt
 
 To reconsile the two pictures let's plot two unit cells in the reciprocal space. (Definition: each of these unit cells is called a **Brillouin zone**, a concept that will come up later.)
 
-![](figures/phonons7.svg)
+```python
+k = np.linspace(0, 2*pi, 300)
+k_dos = np.linspace(0, pi, 20)
+fig, ax = pyplot.subplots()
+ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
+ax.plot(k, dispersion_2m(k), label='acoustic')
+omega_max = dispersion_2m(0, acoustic=False)
+ax.plot(k, omega_max * np.sin(k/4), label='equal masses')
+ax.set_xlabel('$ka$')
+ax.set_ylabel(r'$\omega$')
+ax.set_xticks([0, pi, 2*pi])
+ax.set_xticklabels(['$0$', r'$\pi/2a$', r'$\pi/a$'])
+ax.set_yticks([], [])
+ax.set_xlim(-pi/8, 2*pi+.4)
+ax.set_ylim((0, dispersion_2m(0, acoustic=False) + .2))
+ax.legend(loc='lower right')
+pyplot.vlines([pi, 2*pi], 0, 2.2, linestyles='dashed')
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 Doubling the lattice constant "folds" the band structure on itself, doubling all the bands.
 
@@ -289,7 +433,7 @@ $\rightarrow \rho_{\rm S}(k)=$ (no. of $k$-values / unit of $k$) $=\frac{L}{\pi}
 
 For open boundary conditions $\rightarrow$ running waves $\rightarrow k=$ $...\frac{-4\pi}{L},\frac{-2\pi}{L},0,\frac{2\pi}{L},\frac{4\pi}{L}...$
 
-$\rightarrow \rho_{\rm R}(k)=\frac{L}{2\pi}$, which is lower than for the case of closed boundary conditions, however, in this case the entire $k$-space is filled. 
+$\rightarrow \rho_{\rm R}(k)=\frac{L}{2\pi}$, which is lower than for the case of closed boundary conditions, however, in this case the entire $k$-space is filled.
 
 ## Summary