diff --git a/src/1_einstein_model.md b/src/1_einstein_model.md
index df03afd16b8658a66a11d9c11e6921e18a602913..9f2a291f19d71f28f7cd3491d4402a1afa345ca7 100644
--- a/src/1_einstein_model.md
+++ b/src/1_einstein_model.md
@@ -117,3 +117,52 @@ ax.set_xlabel('$T[K]$')
 ax.set_ylabel('$C/k_B$')
 ax.set_ylim((0, 3));
 ```
+
+## Exercises
+
+### Exercise 1: Heat capacity of a classical oscillator.
+
+Let's refresh the connection of this topic to statistical physics.
+You will need to look up the definition of partition function and how to use it to compute expectation values.
+
+Consider a 1D simple harmonic oscillator with mass $m$ and spring constant $k$.
+The Hamiltonian is given in the usual way by:
+$$
+H = \frac{p^2}{2m}+\frac{k}{2}x^2.
+$$
+
+1. Compute the classical partition function using the following expression:
+$$
+Z = \int_{-\infty}^{\infty}dp \int_{-\infty}^{\infty} dx e^{-\beta H(p,x)}.
+$$
+2. Using the solution of 1., compute the expectation value of the energy, and the expectation value of .
+3. Compute the heat capacity. Check that you get the law of Dulong-Petit but with a different prefactor. 
+4. Explain the difference in the prefactor by considering the number of degrees of freedom.
+
+
+### Exercise 2: Quantum harmonic oscillator
+Consider a 1D quantum harmonic oscillator. Its eigenstates are:
+
+$$
+E_n = \hbar\omega(n+\frac{1}{2}),
+$$
+
+1. Sketch the wave function of this harmonic oscillator for $n=3$.
+2. Compute the quantum partition function using the following expression:
+$$
+Z = \sum_j e^{-\beta E_j}.
+$$
+3. Using the partition function, compute the expectation value of the energy.
+4. Compute the heat capacity. Check that in the high temperature limit you get the same result as in Exercise 1.1.
+  - What temperature can be considered high?
+  - What is the expectation value of $n$?
+
+
+### Exercise 4. Total heat capacity of a diatomic material
+
+Naturally occurring lithium has [two stable isotopes](https://en.wikipedia.org/wiki/Isotopes_of_lithium): $^6$Li  (7.5%) and $^7$Li (92.5%). Let us extend the Einstein model to take into account the different masses of different isotopes.
+
+1. Assume that the strength of the returning force $k$ experienced by each atom is the same. What is the difference in the oscillation frequencies of different isotopes of lithium in the lithium crystal?
+2. Write down the total energy of lithium assuming that all $^6$Li atoms are in $n=2$ vibrational state, and all $^7$Li atoms are in $n=4$ vibrational state.
+3. Write down the total energy of lithium at a temperature $T$ by modifying the Einstein model.
+4. Compute the heat capacity of lithium as a function of $T$.
\ No newline at end of file
diff --git a/src/2_debye_model.md b/src/2_debye_model.md
index abd2b2a00b3c2d19220cf9462d3757f5c7541780..b0157ac3405104ef27c45d9c5b3344490bc50771 100644
--- a/src/2_debye_model.md
+++ b/src/2_debye_model.md
@@ -128,3 +128,64 @@ ax.legend(loc='lower right')
 pyplot.hlines([3], 0, 1.5, linestyles='dashed')
 draw_classic_axes(ax, xlabeloffset=0.3)
 ```
+
+## Exercises
+
+### Exercise 1: Debye model: concepts
+1. Describe the concept of k-space. What momenta are allowed in a 2D system with dimensions $L\times L$?
+2. The probability to find an atom of a 1D solid that originally had a position $x$ at a displacement $\delta x$ is shown on this plot:
+
+```python
+def psi_squared(delta_x, x):
+    return delta_x**2 * np.exp(-delta_x**2) * np.sin(4*np.pi*x)**2
+
+x = np.linspace(0, 1, 200)
+delta_x = np.linspace(-2, 2, 200)
+
+pyplot.imshow(psi_squared(delta_x.reshape((-1, 1)), x.reshape((1, -1))), cmap='gist_heat_r', extent=(0, 3, -1, 1))
+pyplot.ylabel(r'$\delta x$')
+pyplot.xlabel(r'$x$')
+pyplot.xticks((0, 3), ('$0$', '$L$'))
+pyplot.yticks((), ())
+cbar = pyplot.colorbar()
+cbar.set_ticks(())
+cbar.set_label(r'$|\psi^2|$')
+```
+
+Describe how many phonons in which $k$-state this solid has.
+Explain your answer.
+
+??? hint
+
+    There are $n=2$ phonons in the state with $k=4\pi/L$ and $n=2$ phonons in a state with $k=-4\pi/L$.
+
+3. Explain the concept of density of states.
+4. Calculate the density of states $g(\omega)$ for a 3D, 2D and 1D systems with linear dispersion $\omega=vk$.
+
+###  Exercise 2: Debye model in 2D
+
+1. State the assumptions of the Debye theory.
+2. Determine the energy of a two-dimensional solid as a function of $T$ using the Debye approximation (the integral can't be solved analytically).
+3. Calculate the heat capacity in the limit of high $T$ (hint: it goes to a constant).
+4. At low $T$, show that $C_V=KT^{n}$. Find $n$. Express $K$ as a definite integral.
+
+###  Exercise 3: Different phonon modes
+During the lecture we derived the low-temperature heat capacity assuming that all the phonons have the same sound velocity $v$. 
+In reality the longitudinal and transverse modes have different sound velocities (see [Wikipedia](https://en.wikipedia.org/wiki/Sound#Longitudinal_and_transverse_waves) for an illustration of different sound wave types).
+
+Assume that there are two types of excitations:
+
+* One longitudinal mode with $\omega = v_\parallel |k|$
+* Two transverse modes with $\omega = v_\bot |k|$
+
+1. Write down the total energy of phonons in this material *(hint: use the same reasoning as in the [Lithium exercise](1_einstein_model.md#exercise-4-total-heat-capacity-of-a-diatomic-material))*.
+2. Verify that at high $T$ you reproduce the Dulong-Petit law.
+3. Compute the behavior of heat capacity at low $T$.
+
+### Exarcise 4: Anisotropic sound velocities
+Suppose now that the velocity is anisotropic ($v_x \neq v_y \neq v_z$) and $\omega = \sqrt{v_x^2 k_x^2 + v_y^2 k_y^2 + v_z^2 k_z^2}$.
+How does this change the Debye result for the heat capacity?
+
+??? hint
+
+    Write down the total energy as an integral over $k$, then change the integration variables so that the spherical symmetry of the integrand is restored.
\ No newline at end of file