+ c_{\text{V}}=\frac{\partial}{\partial T}\, {\color{red}3}\int_0^{k_D}
+ {\color{blue}\frac{k^2}{2\pi^2}dk} \frac{\hbar {\color{green}ck}}
+ {e^{\beta\hbar {\color{green}ck}}-1}=
+ \frac{\partial}{\partial T}\frac{3\hbar c}{2\pi^2}\int_0^{k_D}
+ \frac{k^3}{e^{\beta\hbar ck}-1}dk=
+ \frac{3\hbar c}{2\pi^2}\int_0^{k_D}
+ \frac{k^3e^{\beta\hbar ck}\beta\hbar ck\frac{1}{T}}
+ {(e^{\beta\hbar ck}-1)^2}dk
+ \]
+ Change of variables: $\beta\hbar ck=x$
+ \[
+ \Rightarrow
+ k=\frac{x}{\beta\hbar c} \quad \textrm{, } \quad
+ dk=\frac{1}{\beta\hbar c} dx
+ \]
+ \[
+ c_{\text{V}}=\frac{3\hbar c}{2\pi^2}\int_0^{\Theta_D/T}
+ \frac{x^3e^xx}{T(\beta\hbar c)^3(e^x-1)^2}\frac{dx}{\beta\hbar c}=
+ \frac{3}{2\pi^2T\beta^4\hbar^3 c^3}\int_0^{\Theta_D/T}
+ \frac{x^4e^x}{(e^x-1)^2}dx
+ \]
+ \[
+ \frac{3}{2\pi^2T\beta^4\hbar^3 c^3}=
+ \frac{3k_{\text{B}}}{2\pi^2\beta^3\hbar^3 c^3}=
+ \frac{3k_{\text{B}}T^33n}{\Theta_{\text{D}}^3}
+ \]
+ \[
+ \Rightarrow
+ c_{\text{V}}=9nk_{\text{B}}\left(\frac{T}{\Theta_{\text{D}}}
+ \right)^3\int_0^{\Theta_D/T}\frac{x^4e^x}{(e^x-1)^2}dx