diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index c0aec168371fc86c10ca91b3379038b94910c74f..6ad51b541b5b667c662faaaba42258d7980ac3cc 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -11,17 +11,6 @@ build lectures:
     expire_in: 1 week
   needs: []
 
-compile mathjax:
-  image: gitlab.kwant-project.org:5005/qt-restricted/registry/mathjax-jsdom
-  script:
-    - find site -name "*.html" | xargs -L1 mathjax-jsdom-mod
-  artifacts:
-    paths:
-      - site
-    expire_in: 1 week
-  needs:
-    - build lectures
-
 .prepare_deploy: &prepare_deploy
   only:
     - branches@solidstate/lectures
@@ -47,7 +36,7 @@ compile mathjax:
   script:
     - "rsync -rv site/* solidstate@tnw-tn1.tudelft.net:$DEPLOY_PATH"
   needs:
-    - compile mathjax
+    - build lectures
 
 deploy delft version:
   <<: *prepare_deploy
diff --git a/docs/10_xray.md b/docs/10_xray.md
index 89f5cd26b112b43af5e435ed39e67632fc71d09c..060feb8d5c10789d415a5f2e9773526dbe6f0ea5 100644
--- a/docs/10_xray.md
+++ b/docs/10_xray.md
@@ -379,11 +379,11 @@ $$
 
 The corresponding phase difference is: 
 
-\begin{align}
+\begin{align*}
 \Delta \phi &= \lvert\mathbf{k} \rvert(\Delta x_1+\Delta x_2)\\
 &= \lvert\mathbf{k}\rvert \lvert\mathbf{R}\rvert(\cos(\theta)+\cos(\theta'))\\
 &= \mathbf{k'}\cdot \mathbf{R} - \mathbf{k}\cdot \mathbf{R} = (\mathbf{k'} - \mathbf{k}) \cdot \mathbf{R}.
-\end{align}
+\end{align*}
 
 However, that is only a phase difference between waves scattered off of two atoms.
 To find the outgoing wave's amplitude, we must sum over the scattered waves from each and every atom in the lattice:
@@ -419,10 +419,10 @@ In order to keep track of the atoms, we define $\mathbf{r}_j$ to be the location
 The distance $\mathbf{r}_j$ is defined with respect to the lattice point from which we construct the unit cell.
 In order to calculate the amplitude of the scattered wave, we must sum not only over all the lattice points but also over the atoms in a single unit cell:
 
-\begin{align}
+\begin{align*}
 A &\propto \sum_\mathbf{R} \sum_j f_j \mathrm{e}^{i\left(\mathbf{G}\cdot(\mathbf{R}+\mathbf{r}_j)-\omega t\right)}\\
 &= \sum_\mathbf{R}\mathrm{e}^{i\left(\mathbf{G}\cdot\mathbf{R}-\omega t\right)}\sum_j f_j\ \mathrm{e}^{i\mathbf{G}\cdot\mathbf{r}_j}
-\end{align}
+\end{align*}
 
 where $f_j$ is the scattering amplitude off of a single atom, and it is called the *form factor*.
 The form factor mainly depends on the chemical element, nature of the scattered wave, and finer details like the electrical environment of the atom.
@@ -455,20 +455,20 @@ As a demonstration of how it happens, let us compute the structure factor of the
 The basis of the conventional FCC unit cell contains four identical atoms.
 With respect to the reference lattice point, these attoms are located at
 
-\begin{align}
+\begin{align*}
 \mathbf{r}_1&=(0,0,0)\\
 \mathbf{r}_2&=\frac{1}{2}(\mathbf{a}_1+\mathbf{a}_2)\\
 \mathbf{r}_3&=\frac{1}{2}(\mathbf{a}_2+\mathbf{a}_3)\\
 \mathbf{r}_4&=\frac{1}{2}(\mathbf{a}_3+\mathbf{a}_1),
-\end{align}
+\end{align*}
 
 with $f_1=f_2=f_3=f_4\equiv f$. Let the reciprocal lattice be described by $\mathbf{G}=h\mathbf{b}_1+k\mathbf{b}_2+l\mathbf{b}_3$, where $h$, $k$ and $l$ are integers. Using the basis described above and the reciprocal lattice, we calculate the structure factor:
 
-\begin{align}
+\begin{align*}
 S&=f\left(\mathrm{e}^{i\mathbf{G}\cdot\mathbf{r}_1}+\mathrm{e}^{i\mathbf{G}\cdot\mathbf{r}_2}+\mathrm{e}^{i\mathbf{G}\cdot\mathbf{r}_3}+\mathrm{e}^{i\mathbf{G}\cdot\mathbf{r}_4}\right)\\
 &=f\left(1+\mathrm{e}^{i(h\mathbf{b}_1\cdot\mathbf{a}_1+k\mathbf{b}_2\cdot\mathbf{a}_2)/2}+\mathrm{e}^{i(k\mathbf{b}_2\cdot\mathbf{a}_2+l\mathbf{b}_3\cdot\mathbf{a}_3)/2}+\mathrm{e}^{i(h\mathbf{b}_1\cdot\mathbf{a}_1+l\mathbf{b}_3\cdot\mathbf{a}_3)/2}\right)\\
 &=f\left(1+\mathrm{e}^{i\pi(h+k)}+\mathrm{e}^{i\pi(k+l)}+\mathrm{e}^{i\pi(h+l)}\right).
-\end{align}
+\end{align*}
 
 Because $h$, $k$ and $l$ are integers, all exponents are either $+1$ or $-1$, and the interference is only present if
 
@@ -539,18 +539,18 @@ plt.axis('off');
 To deduce the values of $\theta$ of a specific crystal, let us put the Laue condition into a more practical form.
 We first take the modulus squared of both sides:
 
-\begin{align}
+\begin{align*}
 \left|\mathbf{G}\right|^2 &= \left|\mathbf{k'}-\mathbf{k} \right|^2 \\
 G^2 &=  2k^2-2\mathbf{k'} \cdot \mathbf{k},
-\end{align}
+\end{align*}
 
 where we used $|\mathbf{k'}| = |\mathbf{k}|$.
 We then substitute the Laue condition $\mathbf{k'} = \mathbf{k}+\mathbf{G}$:
 
-\begin{align}
+\begin{align*}
 \lvert \mathbf{G} \rvert ^2 &= 2k^2-2 \left(\mathbf{k}+\mathbf{G}\right) \cdot \mathbf{k} \\
 &= -2 \mathbf{G} \cdot \mathbf{k}.
-\end{align}
+\end{align*}
 
 Using $\mathbf{k} \cdot \mathbf{G} = \lvert \mathbf{k} \rvert \lvert \mathbf{G} \rvert \cos(\phi)$,  we obtain
 
diff --git a/docs/10_xray_solutions.md b/docs/10_xray_solutions.md
index fa807e8a3ae19c1c5b831a75c28c03351c9f6adc..1325299e6bcdd0ecb13c3ae4d6d25a8d56539410 100644
--- a/docs/10_xray_solutions.md
+++ b/docs/10_xray_solutions.md
@@ -24,21 +24,21 @@ Because we sum over all lattice points, each argument has a different phase. Sum
 
 1. A possible set of BCC primitive lattice vectors is:
 
-    \begin{align}
+    \begin{align*}
     \mathbf{a_1} & = \frac{a}{2} \left(-\hat{\mathbf{x}}+\hat{\mathbf{y}}+\hat{\mathbf{z}} \right) \\
     \mathbf{a_2} & = \frac{a}{2} \left(\hat{\mathbf{x}}-\hat{\mathbf{y}}+\hat{\mathbf{z}} \right) \\
     \mathbf{a_3} & = \frac{a}{2} \left(\hat{\mathbf{x}}+\hat{\mathbf{y}}-\hat{\mathbf{z}} \right).
-    \end{align}
+    \end{align*}
 
     We expect that the volume of the primitive unit cell equals $a^3/2$ because the conventional unit cell has volume $a^3$ and contains two lattice points. We confirm this by calculating the volume of the primitive unit cell using $|\mathbf{a_1}\times\mathbf{a_2}\cdot\mathbf{a_3}| = a^3/2$.
 
 2. The corresponding set of reciprocal lattice vectors
 
-    \begin{align}
+    \begin{align*}
     \mathbf{b_1} & = \frac{2 \pi}{a} \left(\hat{\mathbf{y}}+\hat{\mathbf{z}} \right) \\
     \mathbf{b_2} & = \frac{2 \pi}{a} \left(\hat{\mathbf{x}}+\hat{\mathbf{z}} \right) \\
     \mathbf{b_3} & = \frac{2 \pi}{a} \left(\hat{\mathbf{x}}+\hat{\mathbf{y}} \right),
-    \end{align}
+    \end{align*}
 
 3. The reciprocal lattice forms an FCC lattice. Given the lattice vectors, the volume of the conventional unit cell is $(4\pi/a)^3$.
 
diff --git a/docs/11_nearly_free_electron_model_solutions.md b/docs/11_nearly_free_electron_model_solutions.md
index 63de768faf2e6d811003d7e214519c06e26d162b..47c2e0ea25277c2bce995ece41e4a20d5734a8d0 100644
--- a/docs/11_nearly_free_electron_model_solutions.md
+++ b/docs/11_nearly_free_electron_model_solutions.md
@@ -55,9 +55,9 @@ All $k_n$ that differ by an integer multiple of $2\pi/a$ from $k_0$ have the exa
 \phi_n(x)=\frac{1}{\sqrt{\Omega}} \exp\left[i \left(k_0+\frac{2\pi n}{a}\right)x \right]
 \end{equation*}
 
-\begin{equation}
+\begin{equation*}
 E_n=\frac{\hbar^2}{2m}\left(k_0+\frac{2\pi n}{a}\right)^2
-\end{equation}
+\end{equation*}
 
 ### Question 3.
 ```python
@@ -99,34 +99,34 @@ dispersions(5)
 
 First the kinetic term,
 
-\begin{equation}
+\begin{equation*}
 \left\langle\phi_m|\hat{K}|\phi_n\right\rangle=C_m\frac{\hbar^2k_m^2}{2m}
-\end{equation}
+\end{equation*}
 
 And the potential term,
 
-\begin{align}
+\begin{align*}
 \left\langle\phi_m|V(x)|\psi\right\rangle=\sum_{n=-\infty}^{\infty}C_n\int_0^a V(x) e^{-i\frac{2\pi}{a}(m-n)x}dx
-\end{align}
+\end{align*}
 
 then relabel indices and combine both expressions to find the final answer and expression for $\varepsilon_m$ (which is the free electron dispersion).
 
 ### Question 5.
 From the expression for the energy, it is clear that the difference with respect to the free electron model is given by the Fourier component $V_{m-n}$, describing the coupling between two states $m$ and $n$. The question becomes: when does this term contribute significantly? For that we look at two orthogonal states $\phi_n$ and $\phi_m$, and construct the Hamiltonian in the basis ($\phi_n$,$\phi_m$),  
 
-\begin{equation}
+\begin{equation*}
 \hat{H}=
 \begin{pmatrix}
 \dfrac{\hbar^2 k_n^2}{2m} & V_{n-m}\\
 V_{m-n} & \dfrac{\hbar^2 k_m^2}{2m}
 \end{pmatrix}
-\end{equation}
+\end{equation*}
 
 The eigenvalues of this fellow are 
 
-\begin{equation}
+\begin{equation*}
 E=\frac{\hbar^2(k_n^2+k_m^2)}{4m}\pm\sqrt{\frac{\hbar^4}{16m^2}(k_n^2-k_m^2)^2+|V_{n-m}|^2}.
-\end{equation}
+\end{equation*}
 
 This clearly displays that only if $|k_n|\approx |k_m|$, the band structure will be affected (given that the potential is weak, and therefore small). This nicely demonstrates how an avoided crossing arises.
 
@@ -155,10 +155,10 @@ See the lecture notes.
 ### Question 3.
 We split the Hamiltonian into two parts $H=H_n+H_{\overline{n}}$, where $H_n$ describes a particle in a single delta-function potential well, and $H_\hat{n}$ is the perturbation by the other delta functions:
 
-\begin{align}
+\begin{align*}
 H_n = & \frac{-\hbar^2}{2m}\frac{\partial^2}{\partial x^2} - \lambda\delta(x-na) \\
 H_\overline{n} = & - \lambda \sum_{m\neq n}\delta(x-ma)
-\end{align}
+\end{align*}
 
 such that $H_n|n\rangle = \epsilon_0|n\rangle = -\hbar^2\kappa^2/2m |n\rangle$ with $\kappa=m\lambda/\hbar^2$. We can now calculate 
 
@@ -188,10 +188,10 @@ $$
 
 and
 
-\begin{align}
+\begin{align*}
 \langle n-1|H_\overline{n}|n\rangle = & -\kappa \lambda \sum_{m \neq 0 }\int e^{-\kappa|x-a|} \delta(x-ma)  e^{-\kappa|x|} \\
 =& -\kappa \lambda \sum_{m \neq 0 } e^{-\kappa a|m-1|} e^{-\kappa a |m|} =-\kappa \lambda(e^{ka}+e^{-ka}) \sum_{m=1}^{m=\infty} e^{-2\kappa a m}
-\end{align}
+\end{align*}
 
 In the limit $\kappa a \gg 1$ (i.e., where the distance between the delta functions is large compared to the width of the isolated orbitals), the onsite energy becomes
 
diff --git a/docs/12_band_structures_in_higher_dimensions_solutions.md b/docs/12_band_structures_in_higher_dimensions_solutions.md
index 09800333dc9896afbdcb9c2803924c1db6999012..09c1f367b225d8dc8682d4de35580b3106c16925 100644
--- a/docs/12_band_structures_in_higher_dimensions_solutions.md
+++ b/docs/12_band_structures_in_higher_dimensions_solutions.md
@@ -86,12 +86,12 @@ Four in total: $(\pm\pi/a,\pm\pi/a)$.
 ### Question 3.
 
 Define a basis, e.g.
-\begin{align}
+\begin{align*}
     \left|0\right\rangle &= (\pi/a,\pi/a) \\
     \left|1\right\rangle &= (\pi/a,-\pi/a) \\
     \left|2\right\rangle &= (-\pi/a,-\pi/a) \\
     \left|3\right\rangle &= (-\pi/a,\pi/a)
-\end{align}
+\end{align*}
 The Hamiltonian becomes
 
 $$
diff --git a/docs/13_semiconductors.md b/docs/13_semiconductors.md
index 0a9bc9434fa43100f9f5b3f099c65af85293b950..64c6f3c0edabcd584e4684759b9b25e8384bb0bf 100644
--- a/docs/13_semiconductors.md
+++ b/docs/13_semiconductors.md
@@ -78,10 +78,10 @@ We distinguish three different band filling types: filled, empty and partially f
 Despite being two opposite extreme extreme cases, filled and empty bands are very similar.
 For example, both filled and empty bands carry no electric current:
 
-\begin{align}
+\begin{align*}
 j = 2e \frac{1}{2\pi} \int_{-\pi/a}^{\pi/a} v(k) dk = 2e \frac{1}{2\pi \hbar} \int_{-\pi/a}^{\pi/a} \frac{dE}{dk} \times dk = \\
 2e \frac{1}{2\pi \hbar} [E(\pi/a) - E(-\pi/a)] = 0.
-\end{align}
+\end{align*}
 
 On the other hand, a filled band has an equal number of electrons going forwards and backwards which thus cancel and lead to zero current. 
 Similar results apply to many other physical quantities such as heat capacity and magnetization. 
@@ -187,20 +187,20 @@ draw_classic_axes(ax)
 
 Or in other words
 
-\begin{align}
+\begin{align*}
 E_e &= E_c + \frac{\hbar^2k^2}{2m_e},\\
 E_h &= E_{v,h} + \frac{\hbar^2k^2}{2m_h} = -E_{v} + \frac{\hbar^2k^2}{2m_h}.
-\end{align}
+\end{align*}
 
 Here $E_c$ is the energy of an electron at the bottom of the conduction band and $E_v$ is the energy of an electron at the top of the valence band.
 Observe that because we are describing particles in the valence band as holes, $m_h > 0$ and $E_h > -E_v$.
 
 The corresponding density of states of the two types of particles is
 
-\begin{align}
+\begin{align*}
 g(E) &= (2m_e)^{3/2}\sqrt{E-E_c}/2\pi^2\hbar^3,\\
 g_h(E_h) &= (2m_h)^{3/2}\sqrt{E_h+E_v}/2\pi^2\hbar^3.
-\end{align}
+\end{align*}
 
 ??? question "A photon gives a single electron enough energy to move from the valence band to the conduction band. How many particles does this process create?"
     Two: one electron and one hole.
@@ -248,10 +248,10 @@ Since increasing the Fermi level increases the number of electrons and reduces t
 
 Applying the first two steps of the algorithm:
 
-\begin{gather}
+\begin{gather*}
 n_h = \int_{-E_v}^\infty f_h(E_h) g_h(E_h) dE_h = \int_{-E_v}^\infty\frac{(2m_h)^{3/2}}{2\pi^2\hbar^3}\sqrt{E_h+E_v}\frac{1}{e^{(E_h+E_F)/k_BT}+1}dE_h,\\
 n_e = \int_{E_c}^\infty f(E)g_e(E)dE = \int_{E_c}^\infty\frac{(2m_e)^{3/2}}{2\pi^2\hbar^3} \sqrt{E-E_c}\frac{1}{e^{(E-E_F)/k_BT}+1}dE.
-\end{gather}
+\end{gather*}
 
 Note that whenever calculating the hole dependent quantities, we replace all the relevant physical quantities with their hole equivalents.
 Since the hole energy is opposite $E_h = -E$, we replace the Fermi energy $E_F \to -E_F$ and the bottom of the valence band by $E_v \to -E_v$ in the integration limits.
@@ -319,10 +319,10 @@ Earlier, we deduced that empty and filled bands provide no current.
 We finish the analysis by considering partially filled bands of an intrinsic (pristine) semiconductor.
 To calculate the current, we utilize the Drude model and sum the electron and hole contributions:
 
-\begin{gather}
+\begin{gather*}
 j = -n_e e v_e + n_h e v_h \\
  -m_e v_e /\tau_e = -eE;\quad -m_h v_h /\tau_h = eE.
-\end{gather}
+\end{gather*}
 
 We see that despite opposite velocity signs for electrons and holes, they carry electric current in the same direction:
 
diff --git a/docs/13_semiconductors_solutions.md b/docs/13_semiconductors_solutions.md
index 7302550814f5db0e5c98ab99039005a4259b6708..53d891c456460abc6d9e453f2990eec1adf9b571 100644
--- a/docs/13_semiconductors_solutions.md
+++ b/docs/13_semiconductors_solutions.md
@@ -145,10 +145,10 @@ $$
 
 and when replacing we get two coupled equations:
 
-\begin{align}
+\begin{align*}
 \dot{k_x} &= -\frac{e}{m^*}B k_y\\
 \dot{k_y} &= +\frac{e}{m^*}B k_x
-\end{align}
+\end{align*}
 
 The solution to this equation is circular motion of cyclotron frequency of $\omega_c = \frac{eB}{m^*}$, where the Lorentz
 force is perpendicular to $\nabla_\mathbf{k} E$.
@@ -240,18 +240,18 @@ where $E_h = -E$.
 
 The electron density in the conduction band is given by
 
-\begin{align}
+\begin{align*}
 n_e &= \int_{E_G}^{\infty} f(E)g_c(E)dE\\
 &\approx \int_{E_G}^{\infty} e^{-\beta (E-\mu)} g_c(E) dE\\
 &=\sqrt{\frac{k_B T}{\pi t_{cb}}}\frac{e^{\beta (\mu - E_G)}}{a}.
-\end{align}
+\end{align*}
 
 Analoguous, we find for the hole density
 
-\begin{align}
+\begin{align*}
 n_h &= \int_{0}^{\infty} f(E_h) g_h(E_h) dE_h\\
 &\approx \sqrt{\frac{k_B T}{\pi t_{vb}}}\frac{e^{-\beta \mu}}{a}.
-\end{align}
+\end{align*}
 
 #### Question 5.
 
diff --git a/docs/14_doping_and_devices.md b/docs/14_doping_and_devices.md
index 4d5e789c67b3f09308c9bea21a1d562d826e5c07..b246bb9b26d98bc7af0476076f6a797c2060383b 100644
--- a/docs/14_doping_and_devices.md
+++ b/docs/14_doping_and_devices.md
@@ -182,10 +182,10 @@ When $|N_D-N_A| \gg n_i$ the semiconductor is **extrinsic**. In this case, if $N
 We can now easily find the Fermi level.
 From the first approximation, we know that the simplified relation between $n_{e/h}$ and $E_F$ is:
 
-\begin{align}
+\begin{align*}
 n_e &\approx N_C e^{-(E_c - E_F)/kT},\\
 n_h &\approx N_V e^{(E_v-E_F)/kT}.
-\end{align}
+\end{align*}
 
 We substitute the solutions for $n_e$ and $n_h$ found above and solve for the Fermi level:
 
diff --git a/docs/14_doping_and_devices_solutions.md b/docs/14_doping_and_devices_solutions.md
index 37c4b5429b7ce77202f79fef9d1e326be1c0afac..e96841213709c21983cd31f4034492d94a278921 100644
--- a/docs/14_doping_and_devices_solutions.md
+++ b/docs/14_doping_and_devices_solutions.md
@@ -30,10 +30,10 @@ That means that we are near the transition between extrinsic and intrinsic regim
 In this regime we can neglect $n_D$ and $n_A$ in this exercise, just like in the lecture.
 Writing $n_e n_h = n_i^2$ and $n_e - n_h = N_D - N_A$ and solving these together, we obtain
 
-\begin{align} 
+\begin{align*} 
 n_{e} = \frac{1}{2}(\sqrt{D^2+4n_i^2}+D),\\
 n_{h} = \frac{1}{2}(\sqrt{D^2+4n_i^2}-D)
-\end{align}
+\end{align*}
 
 where $D = N_D - N_A$.
 
@@ -105,27 +105,27 @@ $$ I_s(T) \propto e^{-E_{gap}/k_BT}$$
 ### Question 2.
 This a "particle in a box" problem.
 
-\begin{align}
+\begin{align*}
 -\frac{\hbar^2}{2m_e^{\ast}} \nabla^2 \Psi_e &= (E_e-E_c)\Psi_e\\
 -\frac{\hbar^2}{2m_h^{\ast}} \nabla^2 \Psi_h &= (-E_h-E_v)\Psi_h 
-\end{align}
+\end{align*}
 
 Here and below $E_h$ is the energy of a hole in the valence band.
 
 ### Question 3.
 
-\begin{align}
+\begin{align*}
 E_e = E_c + \frac{\hbar^2}{2m_e^{\ast}} ((\frac{\pi n}{L})^2+k_x^2+k_y^2),\\
 -E_h = E_v - \frac{\hbar^2}{2m_h^{\ast}} ((\frac{\pi n}{L})^2+k_x^2+k_y^2)
-\end{align}
+\end{align*}
 
 ### Question 4.
 This is a 2D electron/hole gas, therefore the DOS per unit area expression is the same as in 2D parabolic dispersion:
 
-\begin{align}
+\begin{align*}
 g_e = \frac{m_e^{\ast}}{\pi\hbar^2},\\
 g_h = \frac{m_h^{\ast}}{\pi\hbar^2}
-\end{align}
+\end{align*}
 
 ### Question 5.
 $L$ can be found here using previous questions, by setting:
diff --git a/docs/1_einstein_model.md b/docs/1_einstein_model.md
index 89b605ebeaa33c81c5f1b98e4769bdcb18fdca02..5b252d17dde2d6bb5d9c57701cdc6319dddc37fc 100644
--- a/docs/1_einstein_model.md
+++ b/docs/1_einstein_model.md
@@ -252,10 +252,10 @@ $$
 
 Substituting $n_B$ into the expression for the oscillator energy, we obtain the expectation value of the energy stored in the oscillator at temperature $T$—the *thermal energy* (which, for brevity, we will simply denote as the total energy $E$):
 
-\begin{align}
+\begin{align*}
 E(T) &= (n_B(\omega_0,T) + \frac{1}{2})\hbar\omega_0 \\
 &= \frac{\hbar\omega_0}{e^{\hbar\omega_0/k_B T}-1} + \frac{1}{2}\hbar\omega_0
-\end{align}
+\end{align*}
 
 A plot of the Bose-Einstein distribution and the expected value of the energy are shown below.
 
@@ -297,12 +297,12 @@ Because the energy in the oscillator becomes approximately constant when $k_B T\
 
 Having found an expression for $E(T)$, we now calculate the heat capacity per harmonic oscillator explicitly by using its definition:
 
-\begin{align}
+\begin{align*}
 C(T) &\equiv \frac{d E(T) }{d T}\\
 &= -\frac{\hbar\omega_0}{\left(e^{\hbar\omega_0/k_B T}-1\right)^2}\frac{d}{d T}\left(e^{\hbar\omega_0/k_B T}-1\right)\\
 &= \frac{\hbar^2\omega_0^2}{k_B T^2}\frac{ e^{\hbar\omega_0/k_B T}}{\left(e^{\hbar\omega_0/k_B T}-1\right)^2}\\
 &=k_B \left(\frac{\hbar\omega_0}{k_B T}\right)^2\frac{ e^{\hbar\omega_0/k_B T}}{\left(e^{\hbar\omega_0/k_B T}-1\right)^2}
-\end{align}
+\end{align*}
 
 We can rewrite this equation into
 
diff --git a/docs/2_debye_model.md b/docs/2_debye_model.md
index c3b03176782ee644170bc2fd87632f106fbabeeb..23aaddf75199f0a122e7959988c920e57ce40957 100644
--- a/docs/2_debye_model.md
+++ b/docs/2_debye_model.md
@@ -122,10 +122,10 @@ where $v_s$ is the _sound velocity_ of a material.
 As discussed in the previous lecture, the quantum mechanical excitations of a harmonic oscillator are called *phonons*, and the expected number of phonons in the oscillator at temperature $T$ is given by the Bose-Einstein distribution $n_B(\beta \hbar \omega(\mathbf{k}))$. Instead of having $3N$ oscillators with the same frequency $\omega_0$ as in the Einstein model, we now have $3N$ oscillators (the vibrational modes) with frequencies depending on $\textbf{k}$ through the dispersion relation $\omega(\mathbf{k}) = v_s|\mathbf{k}|$. Apart from this crucial difference with the Einstein model, the calculation of the expectation value of the energy stored in the vibrational modes proceeds in the same way: 
 The expected value of the total energy stored in the oscillators (which, from now on, we will simply denote as the total energy $E$) is given by the sum of the energy stored in the individual oscillators. These oscillators are characterized by their wavevector $\mathbf{k}$:
 
-\begin{align}
+\begin{align*}
 E &= 3 \sum_\mathbf{k} \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\hbar \omega(\mathbf{k}) n_{B}(\beta \hbar \omega(\mathbf{k}))\right)\\
 &= 3 \sum_\mathbf{k} \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\frac{\hbar\omega(\mathbf{k})}{ e^{\hbar\omega(\mathbf{k})/{k_BT}}-1}\right).
-\end{align}
+\end{align*}
 
 Here we used that the expected occupation number is $n_B(\beta \hbar \omega(\mathbf{k}))$.
 
@@ -195,10 +195,10 @@ $$
 
 We can use this approximation to rewrite the total energy as an integral:
 
-\begin{align}
+\begin{align*}
 E &=  3 \sum_\mathbf{k} \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\frac{\hbar\omega(\mathbf{k})}{ e^{\hbar\omega(\mathbf{k})/k_B T}-1}\right) \\
 &\approx 3 \left(\frac{L}{2\pi}\right)^3 \int \textrm{d} \textbf{k} \left(\frac{1}{2}\hbar\omega(\mathbf{k})+\frac{\hbar\omega(\mathbf{k})}{e^{\hbar\omega(\mathbf{k})/k_B T}-1}\right).
-\end{align}
+\end{align*}
 
 Where $\omega(\mathbf{k}) = v_s |\mathbf{k}|$ is the dispersion relation, and the integral goes over 3-dimensional $k$-space.
 
@@ -213,10 +213,10 @@ $$
 
 Performing the change of variables, we obtain the expression for the total energy in spherical coordinates is
 
-\begin{align}
+\begin{align*}
 E &= 3 \frac{L^3}{(2\pi)^3}\int_0^\infty  4 \pi k^2 \left(\frac{1}{2}\hbar\omega(k)+\frac{\hbar\omega(k)}{ {e}^{\hbar\omega(k)/{k_BT}}-1}\right) \textrm{d}k \\
 &=  3 \frac{L^3}{(2\pi)^3}\frac{4 \pi}{v_s^3}\int_0^\infty \omega^2 \left(\frac{1}{2}\hbar\omega+\frac{\hbar\omega}{ e^{\hbar\omega/{k_BT}}-1}\right) \textrm{d}\omega.
-\end{align}
+\end{align*}
 
 We utilized the dispersion relation $\omega(k) = v_s |k|$ and omitted the absolute value of $k$ due to the integral over $k$ only running from 0 to $\infty$ after conversion to spherical coordinates.
 The integral above can be split up into two factors.
@@ -266,10 +266,10 @@ Despite $E_{\textrm Z}$ diverging towards infinity, does not contribute to $C$.
 The integral depends on the temperature through the $e^{\hbar\omega/k_BT}$ term.
 In order evaluate the integral, we substitute $x\equiv\frac{\hbar\omega}{k_BT}$ and remove the temperature dependence of the integrand:
 
-\begin{align}
+\begin{align*}
 E &= E_Z + \frac{3L^3}{2\pi^2 v_s^3}\frac{\left(k_BT\right)^4}{\hbar^3}\int\limits_0^\infty\frac{x^3}{ {\textrm e}^x-1} {\textrm{d}}x \\
 &= E_Z + \frac{3L^3}{2\pi^2 v_s^3}\frac{\left(k_BT\right)^4}{\hbar^3}\frac{\pi^4}{15},
-\end{align}
+\end{align*}
 
 where we used the fact that the integral is equal to $\frac{\pi^4}{15}$[^4].
 As we can see, the energy scales as $T^4$.
@@ -289,11 +289,11 @@ These modes have $\hbar \omega(\mathbf{k}) \lesssim k_B T$.
     Therefore the number of excited modes is proportional to the volume of a sphere $V_{\textbf{k}} = \frac{4 \pi}{3} |\mathbf{k}|^3$, multiplied by the density of modes in $k$-space, $\left(\frac{L}{2 \pi}\right)^3$. 
     Thus the total number of excited modes is 
 
-    \begin{align}
+    \begin{align*}
     N_\textrm{modes} &= V_{\textbf{k}} \left(\frac{L}{2\pi}\right)^3\\
     &\sim \left( |\mathbf{k}| L \right)^3\\
     &\sim (k_B T L/\hbar v_s)^3
-    \end{align}
+    \end{align*}
 
     where we have substituted $|\mathbf{k}| \simeq k_B T /\hbar v_s$ and left out all numerical factors.
 
@@ -301,10 +301,10 @@ These modes have $\hbar \omega(\mathbf{k}) \lesssim k_B T$.
     Therefore, each mode contributes $k_B$ to the heat capacity (Equipartition theorem). 
     As a result, the heat capacity is
 
-    \begin{align}
+    \begin{align*}
     C &= N_{\textrm{modes}} k_B \\
     &\propto k_B (k_B T L/\hbar v_s)^3.
-    \end{align}
+    \end{align*}
 
 
 ## Debye's interpolation for medium $T$
@@ -338,10 +338,10 @@ $$
 Let us now compute $\omega_D$.
 We know that for a 3D system with $N$ atoms has to have exactly $3N$ phonon modes.
 
-\begin{align}
+\begin{align*}
 3N &= \int_0^{\omega_D} g(\omega)d\omega = \frac{3L^3}{2 \pi^2 v_s^3} \int_0^{\omega_D}\omega^2 {\textrm{d}}\omega\\
 &= \frac{L^3\omega_D^3}{2\pi^2v_s^3},
-\end{align}
+\end{align*}
 
 which gives us
 
diff --git a/docs/2_debye_model_solutions.md b/docs/2_debye_model_solutions.md
index f1cad56b3db6d3fe56a2730d3883bc68dc4b4c87..6c74132a2aa73106541f6ed9a2b354420f6c9496 100644
--- a/docs/2_debye_model_solutions.md
+++ b/docs/2_debye_model_solutions.md
@@ -73,19 +73,19 @@ ax.legend();
 2. The distance between nearest-neighbour points in $\mathbf{k}$-space is $2\pi/L$. The density of $\mathbf{k}$-points in 1, 2, and 3 dimensions is $L/(2\pi)$, $L^2/(2\pi)^2$, and $L^3/(2\pi)^3$ respectively.
 3. Express the number of states between frequencies $0<\omega<\omega_0$ as an integral over k-space. Do so for 1D, 2D and 3D. Do not forget the possible polarizations. We assume that in $d$ dimensions there are $d_p$ polarizations.
     
-    \begin{align}
+    \begin{align*}
     N_\text{states, 1D} & = 1_p \frac{L}{2\pi} \int_0^{k_0} 2 dk \\
     N_\text{states, 2D} & = 2_p \left(\frac{L}{2\pi}\right)^2 \int_0^{k_0} 2\pi k dk \\
     N_\text{states, 3D} & = 3_p  \left(\frac{L}{2\pi}\right)^3 \int_0^{k_0} 4\pi k^2 dk 
-    \end{align}
+    \end{align*}
 
 4.   We use $k=\mathbf{k}| = \omega/v_s$ and $dk = d\omega/v_s$ to get 
 
-    \begin{align}
+    \begin{align*}
     N_\text{states, 1D} & = 1_p \frac{L}{2\pi} \int_0^{\omega_0} 2 \frac{1}{v_s} d\omega  := \int_0^{\omega_0} g_{1D}(\omega) d\omega \\
     N_\text{states, 2D} & = 2_p \left(\frac{L}{2\pi}\right)^2 \int_0^{\omega_0} 2\pi \frac{\omega}{v_s^2} d\omega := \int_0^{\omega_0} g_{2D}(\omega) d\omega \\
     N_\text{states, 3D} & = 3_p  \left(\frac{L}{2\pi}\right)^3 \int_0^{\omega_0} 4\pi \frac{\omega^2}{v_s^3} d\omega := \int_0^{\omega_0} g_{3D}(\omega) d\omega
-    \end{align}
+    \end{align*}
 
     The integral boundaries set the frequency region in which you calculate the density of states.
 
@@ -94,10 +94,10 @@ ax.legend();
 ###  Exercise 2: Debye model in 2D.
 1. The energy stored in the vibrational modes of a two-dimensional Debye solid is:
     
-    \begin{align}
+    \begin{align*}
     E & = \int_0^{\omega_D}(n_B(\omega(\mathbf{k}))+\frac{1}{2})\hbar\omega(\mathbf{k}) d\mathbf{k} \\    
     & = \frac{L^2}{\pi v^2\hbar^2\beta^3}\int_{0}^{\beta\hbar\omega_D}\frac{x^2}{e^{x} - 1}dx + E_0
-    \end{align}
+    \end{align*}
 
   
 2. The high-$T$ limit implies $\beta \rightarrow 0$. Therefore, $n_B \approx k_B T/\hbar\omega$, and the integral becomes particularly illuminating:
@@ -123,10 +123,10 @@ ax.legend();
 ###  Exercise 3: Longitudinal and transverse vibrations with different sound velocities
 1. The key idea is that the total energy in the individual harmonic oscillators (the vibrational modes) is the sum of the energies in the individual oscillators: $E = \int_0^{\omega_D}\frac{\hbar \omega g(\omega)}{e^{\beta\hbar\omega} - 1}d\omega + E_Z$, where $g(\omega) = g_\parallel(\omega) + g_\perp(\omega)$. Using 
     
-    \begin{align}
+    \begin{align*}
     g_\parallel & = 1_p  \frac{L^3}{2\pi^2} \frac{\omega^2}{v_\parallel^3}, \\
     g_\perp     & = 2_p  \frac{L^3}{2\pi^2} \frac{\omega^2}{v_\perp^3},
-    \end{align}
+    \end{align*}
 
     we get 
     
@@ -156,16 +156,16 @@ ax.legend();
 ### Exercise 4: Anisotropic sound velocities.
 In this case, the velocity depends on the direction. Note however that, in contrast with the previous exercise, the polarization does not affect the dispersion of the waves. We get
     
-\begin{align}
+\begin{align*}
 E =& 3_p\left(\frac{L}{2\pi}\right)^3 \int \frac{\hbar\omega(\mathbf{k})}{e^{\beta\hbar\omega(\mathbf{k})}-1} dk_x dk_y dk_z + E_Z \\
 =&   3_p\left(\frac{L}{2\pi}\right)^3\frac{1}{v_x v_y v_z} \int \frac{\hbar\kappa}{e^{\beta\hbar\kappa} - 1}d\kappa_x d\kappa_y d\kappa_z + E_Z,
-\end{align}
+\end{align*}
 
 where we made the substitutions $\kappa_x = k_x v_x,\kappa_y = k_y v_y, \kappa_z = k_z v_z$ so that $\omega = \kappa$. Going to spherical coordinates:
 
-\begin{align}
+\begin{align*}
 E &= 3_p \frac{3_p L^3}{2\pi^2}\frac{1}{v_x v_y v_z}\int_{0}^{\kappa_D} \frac{\hbar\kappa^3}{e^{\beta\hbar\kappa} - 1} d\kappa + E_Z  \\
 & = \frac{3_p L^3}{2\pi^2\hbar^3\beta^4}\frac{1}{v_x v_y v_z}\int_0^{\beta\hbar\kappa_D} \frac{x^3}{e^x - 1}dx + E_Z,
-\end{align}
+\end{align*}
 
 Therefore, $C = \frac{dE}{dT} = \frac{6k_B^4 L^3 T^3}{\pi^2\hbar^3}\frac{1}{v_x v_y v_z}\int_0^{\beta\hbar\kappa_D} \frac{x^3}{e^x - 1}dx$. We see that the result is similar to the one with the linear dispersion, the only difference is the factor $1/(v_x v_y v_z)$ instead of $1/v^3$.
diff --git a/docs/3_drude_model.md b/docs/3_drude_model.md
index 883e3338ff3d2a3a91bd68bf4089e852ef92686a..ce3896836fd25c18c264e3c905e0dfe805a49770 100644
--- a/docs/3_drude_model.md
+++ b/docs/3_drude_model.md
@@ -157,11 +157,11 @@ $$
 
 To find the average velocity, we take a weighted average of these two groups of particles:
 
-\begin{align}
+\begin{align*}
 m\mathbf{v}(t+dt)  &= [m\mathbf{v}(t) + F dt]\left(1 - \frac{dt}{\tau}\right) + 0â‹…\frac{dt}{\tau} \\
   &= [m\mathbf{v}(t) + \mathbf{F} dt] \left(1 - \frac{dt}{\tau}\right) \\
                   &= m\mathbf{v}(t) + dt [\mathbf{F} - m\frac{\mathbf{v(t)}}{\tau}] - \frac{\mathbf{F}}{\tau} dt^2
-\end{align}
+\end{align*}
 
 We now neglect the term proportional to $dt^2$ (it vanishes when $dt → 0$).
 Finally, we recognize that $\left[\mathbf{v}(t+dt) - \mathbf{v}(t)\right]/dt = d\mathbf{v}(t)/dt$, which results in
diff --git a/docs/4_sommerfeld_model.md b/docs/4_sommerfeld_model.md
index 159a1a3606114f648f07e4242bde37653844a0b6..ff9fc016ff425955c830b24a56427d86017eebe6 100644
--- a/docs/4_sommerfeld_model.md
+++ b/docs/4_sommerfeld_model.md
@@ -93,10 +93,10 @@ $$
 Here $\beta = 1/k_{B}T$, $\varepsilon$ is the energy, and $\mu$ is the _chemical potential_. The chemical potential is the characteristic energy above which $n_{F}$ goes to zero. 
 We obtain the number of electrons in the system by summing over the occupied states, as determined the by Fermi-Dirac distribution:
 
-\begin{align}
+\begin{align*}
 N &= 2 \sum_{\mathbf{k}} n_{F}(\beta(\varepsilon-\mu))\\
 &= 2 \left( \frac{L}{2 \pi} \right)^3 \int n_{F}(\beta(\varepsilon(\mathbf{k})-\mu))\mathrm{d} \mathbf{k}.
-\end{align}
+\end{align*}
 
 Just like with phonons, we replaced the discrete sum over $\mathbf{k}$ with a volume integral.
 The factor 2 here accounts for the number of distinct electron states per $\mathbf{k}$ value.
@@ -238,26 +238,26 @@ Let us calculate the density of states for a 3D system.
 Because the free-electron dispersion is isotropic[^2], we use spherical coordinates.
 The total number of states below energy $\varepsilon$ is then
 
-\begin{align}
+\begin{align*}
 N_\text{states} &\overset{\mathrm{3D}}{=}2_s \left(\frac{L}{2\pi}\right)^3\int_{\varepsilon(\mathbf{k}) < \varepsilon}\mathrm{d}{\mathbf{k}}\\
 &=2_s \left(\frac{L}{2\pi}\right)^3 4\pi\int k^2\mathrm{d}k\\
 &=\frac{V}{\pi^2}\int k^2\mathrm{d}k,
-\end{align}
+\end{align*}
 
 We rewrite this expression into an integral over energy using the dispersion relation, substituting $k=\hbar^{-1}\sqrt{2m\varepsilon}$ and $\mathrm{d}k=\hbar^{-1}\sqrt{m/2\varepsilon} d\varepsilon$:
 
-\begin{align}
+\begin{align*}
 N_\text{states} &=\frac{V}{\pi^2}\int\frac{2m \varepsilon}{\hbar^3}\sqrt{\frac{m}{2\varepsilon}}\mathrm{d}\varepsilon\\
 &=\frac{Vm^{3/2}}{\pi^2\hbar^3}\int\sqrt{2\varepsilon}\ \mathrm{d}\varepsilon \\
 & = \int g(\varepsilon) \mathrm{d}\varepsilon \\
-\end{align}
+\end{align*}
 
 We thus find the density of states:
 
-\begin{align}
+\begin{align*}
 g(\varepsilon) &= \frac{ \mathrm{d}N_\text{states}}{ \mathrm{d}\varepsilon}\\
 & =\frac{Vm^{3/2}\sqrt{2\varepsilon}}{\pi^2\hbar^3} \propto\sqrt{\varepsilon}
-\end{align}
+\end{align*}
 
 We observe that the density of states for the 3D parabolic free-electron dispersion is proportional to the square root of energy: 
 
@@ -290,11 +290,11 @@ draw_classic_axes(ax, xlabeloffset=.2)
 ## Relation between the Fermi energy and the number of electrons
 The Fermi energy sets the number of electrons in the system $N$. We can see this by expressing $N$ at $T = 0$ as an integral over the density of states:
 
-\begin{align}
+\begin{align*}
 N &= \int \limits_0^{\infty} n_{F}(\beta(\varepsilon-\mu)) g(\varepsilon) \mathrm{d}\varepsilon\\
 &\overset{\mathrm{T = 0}}{=}\int \limits_0^{\varepsilon_F}g(\varepsilon)\mathrm{d}\varepsilon \\
 &\overset{\mathrm{3D}}{=} \frac{V}{3\pi^2\hbar^3}(2m\varepsilon_F)^{3/2}.
-\end{align}
+\end{align*}
 
 Solving this equation for the Fermi energy yields:
 
@@ -387,19 +387,19 @@ Since the base of the triangle is proportional to $k_BT$ and the height is $\sim
 
 These electrons have gained $k_BT$ of thermal energy, such that the total extra energy is
 
-\begin{align}
+\begin{align*}
 E(T) &= E(T = 0) + N_\mathrm{exc}k_BT\\
 &\approx  E(T = 0) + g(\varepsilon_F)k_B^2T^2.
-\end{align}
+\end{align*}
 
 Therefore, the electron heat capacity $C_e$ is
 
-\begin{align}
+\begin{align*}
 C_e &= \frac{ \mathrm{d}E}{ \mathrm{d}T}\\
  &\approx 2 g(ε_F)k_B^2T\\
  &\overset{\mathrm{3D}}{=} 3 Nk_B\frac{T}{T_F}\\
  &\propto T,
-\end{align}
+\end{align*}
 
 where we used $N=\frac{2}{3}\varepsilon_Fg(\varepsilon_F)$ and the Fermi temperature $T_F \equiv \varepsilon_F/k_B$.
 
diff --git a/docs/4_sommerfeld_model_solutions.md b/docs/4_sommerfeld_model_solutions.md
index 7b5b5103c499b3345453eae146aa03d2a9ff0b06..29daf014e670229c05efb1942fed19f08b991382 100644
--- a/docs/4_sommerfeld_model_solutions.md
+++ b/docs/4_sommerfeld_model_solutions.md
@@ -52,19 +52,19 @@ $$
 
     yielding:
 
-    \begin{align}
+    \begin{align*}
     N_\text{states, 1D} & = 2_s \frac{L}{2\pi} \int_0^{k_0} 2 dk, \\
     N_\text{states, 2D} & = 2_s \left(\frac{L}{2\pi}\right)^2 \int_0^{k_0} 2\pi k dk, \\
     N_\text{states, 3D} & = 2_s  \left(\frac{L}{2\pi}\right)^3 \int_0^{k_0} 4\pi k^2 dk.
-    \end{align}
+    \end{align*}
 
 4. To transform the integral from one over k-space to one over energy, we need the dispersion relation. The integral boundaries are $0<\varepsilon<\varepsilon_0$. Using $k=|\mathbf{k}| = \sqrt{2m\varepsilon/\hbar^2}$ and $dk = \sqrt{m/(2\varepsilon\hbar^2)}d\varepsilon $, we get:
 
-    \begin{align}
+    \begin{align*}
     g_{1D}(\varepsilon) &= 2_s \frac{L}{2\pi} \sqrt{\frac{2m}{\hbar^2\varepsilon}},\\
     g_{2D}(\varepsilon) &= 2_s \frac{L^2}{2\pi}\frac{m}{\hbar^2},\\
     g_{3D}(\varepsilon) &= 2_s \frac{L^3}{2\pi^2}\frac{m}{\hbar^2}\sqrt{\frac{2m\varepsilon}{\hbar^2}} .
-    \end{align}
+    \end{align*}
 
 ### Exercise 2: Applying the free electron model to potassium
 
@@ -129,10 +129,10 @@ $$
 
 3\. Since we have $g(\varepsilon)\propto \varepsilon$ and $E_F=0$, the number of electrons that get thermally excited is given by the area of a triangle with base $k_BT$ and height $g(-k_BT)$:
 
-\begin{align}
+\begin{align*}
 n_{ex} &= \frac{1}{2} g(-k_B T) k_B T\\
 &= \frac{L^2 k_B^2T^2}{\pi c^2}.
-\end{align}
+\end{align*}
 
 From this it follows that the energy difference is given by
 
diff --git a/docs/5_atoms_and_lcao.md b/docs/5_atoms_and_lcao.md
index b81e1ec251aa99076d3d3dad578259f76ce24331..233c7795cfcfbdf043e59253f5362425fd8ffc1d 100644
--- a/docs/5_atoms_and_lcao.md
+++ b/docs/5_atoms_and_lcao.md
@@ -140,10 +140,10 @@ Additionally, we assume that the atoms are sufficiently far apart, such that the
 
 If the atoms are far apart from each other such that they do not interact, the eigenstates of the electrons are the atomic orbitals. If we call the atomic orbital of the electron on the first atom $|1\rangle$ and that of the electron on the second atom $|2\rangle$, we have:
 
-\begin{align}
+\begin{align*}
 (\hat{V}_1 + \hat{K})|1\rangle = \varepsilon_0|1\rangle,\\
 (\hat{V}_2 + \hat{K})|2\rangle = \varepsilon_0|2\rangle.
-\end{align}
+\end{align*}
 
 Key idea: to find the wavefunction of the electron on the molecule - the _molecular orbital_ -, we search for a solution that is a *linear combination of the atomic orbitals* (LCAO):
 
@@ -364,10 +364,10 @@ $$|\psi\rangle = \phi_1|1\rangle + \phi_2|2\rangle.$$
 
 The orbitals $|1\rangle$ and $|2\rangle$ correspond to the wave functions of the electron when there is only a single delta peak present:
 
-\begin{align}
+\begin{align*}
 H_1 |1\rangle &= \epsilon_0 |1\rangle,\\
 H_2 |2\rangle &= \epsilon_0 |2\rangle.
-\end{align}
+\end{align*}
 
 We start by calculating the wavefunction of an electron bound to a single delta-function potential welll.
 To do so, you first need to set up the Schrödinger equation for a single electron bound to a delta-function potential well.
diff --git a/docs/5_atoms_and_lcao_solutions.md b/docs/5_atoms_and_lcao_solutions.md
index b70d7a84cc68fac808ab498369dbeb482ccc11ef..ff1a09a2e7610a1df0107177e7fae227580d0074 100644
--- a/docs/5_atoms_and_lcao_solutions.md
+++ b/docs/5_atoms_and_lcao_solutions.md
@@ -35,12 +35,12 @@ search:
 
 3. Based on Madelung's rule, we would expect
 
-    \begin{align}
+    \begin{align*}
     \textrm{Cu} &= [\textrm{Ar}]4s^23d^9\\
     \textrm{Pd} &= [\textrm{Kr}]5s^24d^8\\
     \textrm{Ag} &= [\textrm{Kr}]5s^24d^9\\
     \textrm{Au} &= [\textrm{Xe}]6s^24f^{14}5d^9
-    \end{align}
+    \end{align*}
 
     Deviations from this rule are predominantly caused by the interactions between the electrons, which render some orbitals more energetically favorable and others less so.
 
@@ -94,17 +94,17 @@ search:
 
     The onsite energy is
 
-    \begin{align}
+    \begin{align*}
         E_0 = \langle 1 | \hat{H} | 1\rangle & =  \langle 1 | \hat{V}_1 + \hat{V}_2 + \hat{K}| 1\rangle = E + \langle 1 | \hat{V}_2 | 1\rangle = E -V_0 \int \psi_1(x) \delta(x-x_2) \psi_1(x) dx \\
         =& \varepsilon_0 - V_0 \kappa e^{-2\kappa|x_2-x_1|}
-    \end{align}  
+    \end{align*}  
 
     The hopping is
 
-    \begin{align}
+    \begin{align*}
         -t =  \langle 1 | \hat{H} | 2\rangle & = \langle 1 | \hat{V}_1 + \hat{V}_2 + \hat{K}| 2\rangle =  \langle 1 | \hat{V}_1 | 2\rangle = -V_0 \int \psi_1(x) \delta(x-x_1) \psi_2(x) dx \\
         &=  -V_0 \kappa e^{-\kappa|x_2-x_1|}
-    \end{align}  
+    \end{align*}  
 
     So we get
 
diff --git a/docs/6_bonds_and_spectra.md b/docs/6_bonds_and_spectra.md
index dbf1719e6a00e6f33770b5bd38d0480b469e8f18..a15c5406a3adc85c10e586a6f08321a422f34ad5 100644
--- a/docs/6_bonds_and_spectra.md
+++ b/docs/6_bonds_and_spectra.md
@@ -436,11 +436,11 @@ make_plot(-y_max, deviation_arr, max_m);
 Let us denote the deviation of atom $i$ from its equilibrium position by $u_i$.
 Newton's equations of motion for this system are then given by
 
-\begin{align}
+\begin{align*}
 m \ddot{u}_1 &= - \kappa (u_1 - u_2), \\
 m \ddot{u}_2 &= - \kappa (u_2 - u_1) - \kappa (u_2 - u_3), \\
 m \ddot{u}_3 &= - \kappa (u_3 - u_2).
-\end{align}
+\end{align*}
 
 We write this system of equations in matrix form
 
@@ -493,11 +493,11 @@ Furthermore, we assume hopping only between the nearest neighbors and assume tha
 We also assume that the orbitals are orthogonal to eachother.
 Just as we did in the previous lecture, we use the Schrödinger equation $H |\psi\rangle = E |\psi\rangle$ to set up a system of equations:
 
-\begin{align}
+\begin{align*}
 E \varphi_1 &= E_0 \varphi_1 - t \varphi_2\\
 E \varphi_2 &= E_0 \varphi_2 - t \varphi_1 - t \varphi_3\\
 E \varphi_3 &= E_0 \varphi_3 -t \varphi_2.
-\end{align}
+\end{align*}
 
 Again, we write this in a matrix form:
 
diff --git a/docs/7_tight_binding.md b/docs/7_tight_binding.md
index cb081209d8552ff3d35f535e71264fc276c61a52..36179cad8f8d941a6d178895f4742104aef7bfa1 100644
--- a/docs/7_tight_binding.md
+++ b/docs/7_tight_binding.md
@@ -297,10 +297,10 @@ $$ H = E(k) + e\mathcal{E}x.$$
 
 To derive how particles with an arbitrary dispersion relation move, we recall the [Hamilton's equations](https://en.wikipedia.org/wiki/Hamiltonian_mechanics) for particle velocity $v$ and force $F$:
 
-\begin{align}
+\begin{align*}
 v \equiv \frac{dr}{dt} &= \frac{\partial H(p, r)}{\partial p}\\
 F \equiv \frac{dp}{dt} &= -\frac{\partial H(p, r)}{\partial r}.
-\end{align}
+\end{align*}
 
 Substituting $p = \hbar k$ into the first equation we arrive to the expression for the electron **group velocity** $v \equiv \hbar^{-1}\partial E/\partial k$.
 From the second equation we obtain that the force acting on electron in a band stays $-e\mathcal{E}$, which in turn gives results in the acceleration
diff --git a/docs/8_many_atoms.md b/docs/8_many_atoms.md
index df5c020f05da18eec63306400fbe10443b062197..4d15e5c821ae4380a01fe35a6fb0624b5ed4c79a 100644
--- a/docs/8_many_atoms.md
+++ b/docs/8_many_atoms.md
@@ -52,10 +52,10 @@ In our choice of unit cell, atom 1 has mass $m_1$ and atom 2 has mass $m_2$.
 
 Having specified the degrees of freedom, let's write down the equations of motion:
 
-\begin{align}
+\begin{align*}
 m_1\ddot{u}_{1,n}&=\kappa(u_{2,n}-2u_{1,n}+u_{2,n-1})\\
 m_2\ddot{u}_{2,n}&=\kappa(u_{1, n} - 2u_{2,n}+u_{1,n+1}).
-\end{align}
+\end{align*}
 
 The new ansatz is conceptually the same as before: all unit cells should behave the same up to a phase factor:
 
diff --git a/docs/8_many_atoms_solutions.md b/docs/8_many_atoms_solutions.md
index 975e5bd8b43660a0013b84516bb5460540b97d29..9ff34cfcc1a41c562d032947cb3fcb85c7920315 100644
--- a/docs/8_many_atoms_solutions.md
+++ b/docs/8_many_atoms_solutions.md
@@ -218,11 +218,11 @@ If we have one electron per atom, all states below $E_F = \varepsilon_0$ are fil
 
 2. The equations of motion are
 
-    \begin{align}
+    \begin{align*}
     m\ddot{u}_{n,1} &= -\kappa_1(u_{n,1} - u_{n,2}) - \kappa_3(u_{n,1} - u_{n-1,3})\\
     m\ddot{u}_{n,2} &= -\kappa_2(u_{n,2} - u_{n,3}) - \kappa_1(u_{n,2}-u_{n,1})\\
     m\ddot{u}_{n,3} &= -\kappa_3(u_{n,3} - u_{n+1,1}) - \kappa_2(u_{n,3}-u_{n,2})
-    \end{align}
+    \end{align*}
 
 3. Substitute the following Ansatz into the equations of motion:
 
diff --git a/docs/9_crystal_structure.md b/docs/9_crystal_structure.md
index afb0c7c6221b795fd5d735c1cc6751058c0c5f37..2e2a681ae8de798bcd77faeb3c0fd13a0fb27be2 100644
--- a/docs/9_crystal_structure.md
+++ b/docs/9_crystal_structure.md
@@ -418,10 +418,10 @@ In order to map out the entire lattice, we need to include an extra star in the
 Since there is one star at the corner of the unit cell and one in the centre, the basis is: $\star = (0,0),(1/2,1/2)$.
 The corresponding locations of the stars with respect to the reference lattice point are then given by:
 
-\begin{align}
+\begin{align*}
 \mathbf{R}_{\left[n_{1}, n_{2}\right]}^{\mathrm{corner}} &= n_{1} \mathbf{a}_1 +n_2 \mathbf{a}_2\\
 \mathbf{R}_{\left[n_{1}, n_{2}\right]}^{\mathrm{centre}} &= 1/2\left( \mathbf{a}_1 + \mathbf{a}_2 \right) + n_{1} \mathbf{a}_1 +n_2 \mathbf{a}_2,
-\end{align}
+\end{align*}
 
 for $ n_{1}, n_{2} \in \mathbb{Z}$.
 Alternatively, one can say that the complete crystal structure is made up from two interpenetrating orthogonal lattices - one centred at $[0,0]$ and another at $[1/2,1/2]$, each with basis $\star = (0,0)$.
@@ -631,9 +631,9 @@ $$
 
 We first look at the $y$-component of $\mathbf{r}_1$, which is $a$. Because only $\mathbf{a}_2$ has a nonzero $y$ component, we conclude $f_2=2/3$. We then look at the $x$-component of $\mathbf{r}_1$, which is $\sqrt{3}a$. It should satisfy
 
-\begin{align}
+\begin{align*}
 \sqrt{3}a &= f_1 a_{1,x} + f_2 a_{2,x}\\
-\end{align}
+\end{align*}
 
 which we solve to find $f_1 = 2/3$. Hence, the fractional coordinates of the second atom are $\mathrm{C}(2/3, 2/3)$.
 
@@ -1074,11 +1074,11 @@ $$
 
 We showed earlier that the number of atoms in the fcc unit cell is 4. Therefore, we now have all the information to calculate the filling factor:
 
-\begin{align}
+\begin{align*}
 F &= \frac{ N_{\mathrm{atom}} V_{\mathrm{atom}} }{V_{\mathrm{cell}}}\\
 &= \frac{4 \times \frac{\pi a^3}{12\sqrt{2}}}{a^3}\\
 &= \frac{\pi}{3\sqrt{2}} \approx 0.74
-\end{align}
+\end{align*}
 
 This means that appoximately $74 \%$ of the fcc unit cell is occupied by atoms.
 One might ask if this is the best we can achieve.
diff --git a/docs/9_crystal_structure_solutions.md b/docs/9_crystal_structure_solutions.md
index 847847750833e063c88f8cc158b48daf380a5111..28490bde14064054aa8f7db49f7907abfd491b57 100644
--- a/docs/9_crystal_structure_solutions.md
+++ b/docs/9_crystal_structure_solutions.md
@@ -16,22 +16,22 @@ from math import pi
 
 3. Let us first consider the FCC lattice. Its primitive lattice vectors vectors are
 
-    \begin{align}
+    \begin{align*}
     \mathbf{a}_1 &= \frac{a}{2}(\mathbf{\hat{x}} + \mathbf{\hat{y}})\\
     \mathbf{a}_2 &= \frac{a}{2}(\mathbf{\hat{x}} + \mathbf{\hat{z}})\\
     \mathbf{a}_3 &= \frac{a}{2}(\mathbf{\hat{y}} + \mathbf{\hat{z}}).
-    \end{align}
+    \end{align*}
 
     With respect to the conventional unit cell, the basis in fractional coordinates is $\bigcirc(1/2,1/2,0)$, $\bigcirc(1/2,0,1/2)$, $\bigcirc(0,1/2,1/2)$ and $\bigcirc(0,0,0)$.
     With respect to the primitive unit cell, the basis is $\bigcirc(0,0,0)$.
     Let us now consider the BCC lattice.
     The primitive lattice vectors are
 
-    \begin{align}
+    \begin{align*}
     \mathbf{a}_1 &= \frac{a}{2}(\mathbf{\hat{x}} + \mathbf{\hat{y}} + \mathbf{\hat{z}})\\
     \mathbf{a}_2 &= a\mathbf{\hat{x}}\\
     \mathbf{a}_3 &= a\mathbf{\hat{y}}.
-    \end{align}
+    \end{align*}
 
     The basis of the conventional unit cell is $\bigcirc(0,0,0)$ and $\bigcirc(1/2,1/2,1/2)$.
     For the primitive unit cell the basis is $\bigcirc(0,0,0)$.
@@ -106,11 +106,11 @@ If the filled and empty circles are identical particles, the nearest-neighbour d
 
 1. The conventional unit cell of diamond consists of two intertwined fcc lattices. A possible set of primitive lattice vectors is
 
-    \begin{align}
+    \begin{align*}
     \mathbf{a_1} &= \frac{a}{2} \left(\hat{\mathbf{x}}+\hat{\mathbf{y}} \right) \\
     \mathbf{a_2} &= \frac{a}{2} \left(\hat{\mathbf{x}}+\hat{\mathbf{z}} \right) \\
     \mathbf{a_3} &= \frac{a}{2} \left(\hat{\mathbf{y}}+\hat{\mathbf{z}} \right).
-    \end{align}
+    \end{align*}
 
     The volume of the primitive unit cell is
 
diff --git a/docs/extra_exercises_solutions.md b/docs/extra_exercises_solutions.md
index 62367073f217df2c1cb5826c7bf6e8f68169a89f..991c9fa56d331193c1827f7210f260b2a404769a 100644
--- a/docs/extra_exercises_solutions.md
+++ b/docs/extra_exercises_solutions.md
@@ -12,27 +12,27 @@
 
     The factor 2 is due to positive and negative $k$-values having equal enery
 
-    \begin{align}
+    \begin{align*}
     g_{2D}(k)\textrm{d} k &= \left(\frac{L}{2\pi}\right)^2 2\pi k \mathrm{d} k,\\
     g_{3D}(k)\textrm{d} k &= \left(\frac{L}{2\pi}\right)^3 4\pi k^2 \mathrm{d} k
-    \end{align}
+    \end{align*}
 
 3. 
 
-    \begin{align}
+    \begin{align*}
     g(k)\textrm{d} k &= \left(\frac{L}{2\pi}\right)^n S_{n-1}(k)\textrm{d} k \\
     &= \left(\frac{L}{2\pi}\right)^n \frac{2\pi^{\frac{n}{2}}k^{n-1}}{\Gamma(\frac{n}{2})}\textrm{d} k
-    \end{align}
+    \end{align*}
 
 4. See hint of the question
 
 5. 
 
-    \begin{gather}
+    \begin{gather*}
     g(\varepsilon)\textrm{d} \varepsilon = 2_s g(k)  \textrm{d} k\\
     g(\varepsilon)=2_sg(k(\varepsilon))\frac{\textrm{d} k}{\textrm{d} \varepsilon}\\
     g(\varepsilon) = \frac{2_s}{\Gamma(\frac{n}{2})}\left(\frac{L}{\hbar}\sqrt{\frac{m}{2\pi}}\right)^n (\varepsilon)^{\frac{n}{2}-1}
-    \end{gather}
+    \end{gather*}
 
 6. 
 
@@ -60,10 +60,10 @@ $$
 
 3.
 
-\begin{align}
+\begin{align*}
 E(T)-E(T=0) &= \frac{\pi^2}{6}(k_B T)^2\frac{\partial}{\partial \varepsilon}\left(\varepsilon g(\varepsilon)\right)\bigg|_{\varepsilon=\varepsilon _F}\\
 &\approx 8.356 10^8 eV
-\end{align}
+\end{align*}
 
 5.
 $C_v = 1.6713.10^6 eV/K$
diff --git a/docs/scripts/katex.js b/docs/scripts/katex.js
new file mode 100644
index 0000000000000000000000000000000000000000..70683b41defdd7da65fb3906443cb579553bdad9
--- /dev/null
+++ b/docs/scripts/katex.js
@@ -0,0 +1,11 @@
+document$.subscribe(({ body }) => {
+    renderMathInElement(body, {
+      delimiters: [
+        { left: "$$",  right: "$$",  display: true },
+        { left: "$",   right: "$",   display: false },
+        { left: "\\(", right: "\\)", display: false },
+        { left: "\\[", right: "\\]", display: true }
+      ],
+    })
+  })
+
diff --git a/docs/scripts/mathjaxhelper.js b/docs/scripts/plotlyhelper.js
similarity index 100%
rename from docs/scripts/mathjaxhelper.js
rename to docs/scripts/plotlyhelper.js
diff --git a/mkdocs.yml b/mkdocs.yml
index 9f34cc4398bbb02b28fb7299bdb7c49b9217099c..7d0b9a8bb0d5f58e2d45380ba35dd45955c6cfec 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -62,10 +62,6 @@ theme:
 plugins:
   - search
   - social
-  - mkdocs-simple-hooks:
-      hooks:
-        on_page_read_source: "execute:on_page_read_source"
-        on_page_markdown: "execute:on_page_markdown"
 
 markdown_extensions:
   - pymdownx.arithmatex:
@@ -81,15 +77,21 @@ markdown_extensions:
   - footnotes
   - meta
 
-extra_css:
-  - 'https://use.fontawesome.com/releases/v5.8.1/css/all.css'
-  - 'styles/mathjax_hotfix.css'
-  - 'styles/thebelab.css'
+hooks:
+  - execute.py
 
 extra_javascript:
-  - 'scripts/mathjaxhelper.js'
+  - scripts/katex.js
+  - scripts/plotlyhelper.js
+  - https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.7/katex.min.js
+  - https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.7/contrib/auto-render.min.js
   - 'https://polyfill.io/v3/polyfill.min.js?features=es6'
   - 'https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js'
   - 'https://cdn.plot.ly/plotly-latest.min.js'
 
+extra_css:
+  - https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.16.7/katex.min.css
+  - 'https://use.fontawesome.com/releases/v5.8.1/css/all.css'
+  - 'styles/thebelab.css'
+
 copyright: "Copyright © 2017-2022 Delft University of Technology, CC-BY-SA 4.0."
diff --git a/requirements.txt b/requirements.txt
index 427e333d88592ca631acfd3e8bde3a6da3d0dde2..cdafd21059fadff47364a5d1968901942a86f911 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,6 +2,5 @@ wikitables
 pandas
 python-markdown-math
 mkdocs-material
-mkdocs-simple-hooks
 pillow
 cairosvg