Peeter Joot's (OLD) Blog.

• Recent Comments

 Adam C Scott on avoiding gdb signal noise… Ken on Scotiabank iTrade RESP …… Alan Ball on Oops. Fixing a drill hole in P… Peeter Joot's B… on Stokes theorem in Geometric… Exploring Stokes The… on Stokes theorem in Geometric…

• 301,175

Posts Tagged ‘pressure’

A final pre-exam update of my notes compilation for ‘PHY452H1S Basic Statistical Mechanics’, Taught by Prof. Arun Paramekanti

Posted by peeterjoot on April 22, 2013

Here’s my third update of my notes compilation for this course, including all of the following:

April 21, 2013 Fermi function expansion for thermodynamic quantities

April 20, 2013 Relativistic Fermi Gas

April 10, 2013 Non integral binomial coefficient

April 10, 2013 energy distribution around mean energy

April 09, 2013 Velocity volume element to momentum volume element

April 04, 2013 Phonon modes

April 03, 2013 BEC and phonons

April 03, 2013 Max entropy, fugacity, and Fermi gas

April 02, 2013 Bosons

April 02, 2013 Relativisitic density of states

March 28, 2013 Bosons

plus everything detailed in the description of my previous update and before.

Advertisements

Fermi-Dirac function expansion for thermodynamic quantities

Posted by peeterjoot on April 21, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

In section 8.1 of [1] are some Fermi-Dirac \index{Fermi-Dirac function} expansions for $P$, $N$, and $U$. Let’s work through these in detail.

Our starting point is the relations

\begin{aligned}P V \beta = \ln Z_{\mathrm{G}} = \sum \ln \left( 1 + z e^{-\beta \epsilon} \right)\end{aligned} \hspace{\stretch{1}}(1.0.1a)

\begin{aligned}N = \sum \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }}.\end{aligned} \hspace{\stretch{1}}(1.0.1b)

Recap. Density of states

We’ll employ the 3D non-relativisitic density of states

\begin{aligned}\mathcal{D}(\epsilon) &= \sum_\mathbf{k} \delta(\epsilon - \epsilon_\mathbf{k}) \\ &\sim V \int \frac{d^3 \mathbf{k}}{(2 \pi)^3}\delta(\epsilon - \epsilon_\mathbf{k}) \\ &= \frac{4 \pi V}{(2 \pi)^3}\int dk k^2 \delta\left( \epsilon - \frac{\hbar^2 k^2}{2 m} \right) \\ &= \frac{4 \pi V}{(2 \pi)^3}\int dk k^2 \frac{ \delta\left( k - \sqrt{2 m \epsilon}/\hbar \right)}{ \frac{\hbar^2}{m} \frac{\sqrt{2 m \epsilon}}{\hbar}} \\ &= \frac{2 V}{(2 \pi)^2 }\frac{m}{\hbar^2}\sqrt{\frac{2 m \epsilon}{\hbar^2}},\end{aligned} \hspace{\stretch{1}}(1.0.1b)

or

\begin{aligned}\boxed{\mathcal{D}(\epsilon)=\frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\epsilon^{1/2}.}\end{aligned} \hspace{\stretch{1}}(1.0.1b)

Density

Now let’s make our integral approximation of the sum for $N$. That is

\begin{aligned}N &= g \int d\epsilon \mathcal{D}(\epsilon) \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }} \\ &= g \frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\int_0^\infty d\epsilon \frac{\epsilon^{1/2}}{ z^{-1} e^{\beta \epsilon} + 1 } \\ &= g \frac{V}{(2 \pi)^2 \beta^{3/2}}\left( \frac{2 m}{\hbar^2} \right)^{3/2}\int_0^\infty du \frac{u^{1/2}}{ z^{-1} e^{u} + 1 } \\ &= g \frac{V}{(2 \pi)^2 \beta^{3/2}}\left( \frac{2 m}{\hbar^2} \right)^{3/2}\Gamma(3/2) f_{3/2}(z) \\ &= g \frac{V}{(2 \pi)^2 \beta^{3/2}}\frac{\left( 2 m k_{\mathrm{B}} T \right)^{3/2}}{\hbar^3}\frac{1}{{2}} \sqrt{\pi}f_{3/2}(z)\\ &= g V \not{{2}} \pi\frac{\left( 2 m k_{\mathrm{B}} T \right)^{3/2}}{h^3}\frac{1}{{\not{{2}}}} \sqrt{\pi}f_{3/2}(z),\end{aligned} \hspace{\stretch{1}}(1.0.1b)

or

\begin{aligned}\frac{N}{V} = g \frac{\left( 2 \pi m k_{\mathrm{B}} T \right)^{3/2}}{h^3}f_{3/2}(z).\end{aligned} \hspace{\stretch{1}}(1.0.5)

With

\begin{aligned}\lambda = \frac{h}{\sqrt{ 2 \pi m k_{\mathrm{B}} T }},\end{aligned} \hspace{\stretch{1}}(1.0.6)

this gives us the desired density result from the text

\begin{aligned}\boxed{\frac{N}{V}=\frac{g}{\lambda^3} f_{3/2}(z).}\end{aligned} \hspace{\stretch{1}}(1.0.7)

Pressure

For the pressure, we can do the same, but have to integrate by parts

\begin{aligned}P V \beta &= g \sum \ln \left( 1 + z e^{-\beta \epsilon} \right) \\ &\sim g \frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\int_0^\infty d\epsilon \epsilon^{1/2} \ln \left( 1 + z e^{-\beta \epsilon} \right) \\ &= - g \frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\int_0^\infty d\epsilon \frac{2}{3} \epsilon^{3/2} \frac{-\beta z e^{-\beta \epsilon} }{ 1 + z e^{-\beta \epsilon} } \\ &= g\frac{V}{(2 \pi)^2 }\left( \frac{2 m}{\hbar^2} \right)^{3/2}\frac{2}{3} \frac{1}{{\beta^{3/2}}}\int_0^\infty dx\frac{x^{3/2}}{z^{-1} e^{x} + 1 } \\ &= g\frac{2}{3} 2 \pi V\frac{\left( 2 m k_{\mathrm{B}} T \right)^{3/2}}{h^3 }\Gamma(5/2)f_{5/2}(z) \\ &= g\frac{2}{3} 2 \pi V\frac{\left( 2 m k_{\mathrm{B}} T \right)^{3/2}}{h^3 }\frac{3}{2} \frac{1}{2} \sqrt{\pi}f_{5/2}(z) \\ &= g V\frac{\left( 2 \pi m k_{\mathrm{B}} T \right)^{3/2}}{h^3 }f_{5/2}(z),\end{aligned} \hspace{\stretch{1}}(1.0.7)

or

\begin{aligned}\boxed{P \beta = \frac{g}{\lambda^3} f_{5/2}(z).}\end{aligned} \hspace{\stretch{1}}(1.0.9)

Energy

The average energy is the last thermodynamic quantity to come very easily. We have

\begin{aligned}U &= - \frac{\partial {}}{\partial {\beta}} \ln Z_{\mathrm{G}} \\ &= - \frac{\partial {T}}{\partial {\beta}} \frac{\partial {}}{\partial {T}} \ln Z_{\mathrm{G}} \\ &= - \frac{\partial {(1/k_{\mathrm{B}} T)}}{\partial {\beta}} \frac{\partial {}}{\partial {T}} P V \beta \\ &= \frac{1}{{k_{\mathrm{B}} \beta^2}}\frac{\partial {}}{\partial {T}} \frac{g V}{\lambda^3} f_{5/2}(z) \\ &= g V k_{\mathrm{B}} T^2f_{5/2}(z)\frac{\partial {}}{\partial {T}} \frac{\left( 2 \pi m k_{\mathrm{B}} T \right)^{3/2}}{h^3} \\ &= \frac{3}{2} \frac{g V k_{\mathrm{B}} T}{\lambda^3}f_{5/2}(z).\end{aligned} \hspace{\stretch{1}}(1.0.9)

From eq. 1.0.7, we have

\begin{aligned}\frac{g V}{\lambda^3} = \frac{N}{f_{3/2}(z) },\end{aligned} \hspace{\stretch{1}}(1.0.11)

so the energy takes the form

\begin{aligned}\boxed{U = \frac{3}{2} N k_{\mathrm{B}} T \frac{f_{5/2}(z)}{f_{3/2}(z) }.}\end{aligned} \hspace{\stretch{1}}(1.0.11)

We can compare this to the ratio of pressure to density

\begin{aligned}\frac{P \beta}{n} = \frac{f_{5/2}(z)}{f_{3/2}(z) },\end{aligned} \hspace{\stretch{1}}(1.0.11)

to find

\begin{aligned}U= \frac{3}{2} N k_{\mathrm{B}} T \frac{P V \beta}{N}= \frac{3}{2} P V,\end{aligned} \hspace{\stretch{1}}(1.0.11)

or

\begin{aligned}\boxed{P V = \frac{2}{3} U.}\end{aligned} \hspace{\stretch{1}}(1.0.11)

References

[1] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996.

Relativisitic Fermi gas

Posted by peeterjoot on April 20, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Question: Relativisitic Fermi gas ([1], pr 9.3)

Consider a relativisitic gas of $N$ particles of spin $1/2$ obeying Fermi statistics, enclosed in volume $V$, at absolute zero. The energy-momentum relation is

\begin{aligned}\epsilon = \sqrt{(p c)^2 + \epsilon_0^2},\end{aligned} \hspace{\stretch{1}}(1.1)

where $\epsilon_0 = m c^2$, and $m$ is the rest mass.

Find the Fermi energy at density $n$.

With the pressure $P$ defined as the average force per unit area exerted on a perfectly-reflecting wall of the container.

Set up expressions for this in the form of an integral.

Define the internal energy $U$ as the average $\epsilon - \epsilon_0$.
Set up expressions for this in the form of an integral.

Show that $P V = 2 U/3$ at low densities, and $P V = U/3$ at high densities. State the criteria for low and high densities.

There may exist a gas of neutrinos (and/or antineutrinos) in the cosmos. (Neutrinos are massless Fermions of spin $1/2$.) Calculate the Fermi energy (in eV) of such a gas, assuming a density of one particle per $\text{cm}^3$.

Attempt exact evaluation of the various integrals.

Answer

We’ve found [3] that the density of states associated with a 3D relativisitic system is

\begin{aligned}\mathcal{D}(\epsilon) = \frac{4 \pi V}{(c h)^3} \epsilon \sqrt{\epsilon^2 -\epsilon_0^2},\end{aligned} \hspace{\stretch{1}}(1.0.2)

For a given density $n$, we can find the Fermi energy in the same way as we did for the non-relativisitic energies, with the exception that we have to integrate from a lowest energy of $\epsilon_0$ instead of $0$ (the energy at $\mathbf{p} = 0$). That is

\begin{aligned}n &= \frac{N}{V} \\ &= \left( 2 \frac{1}{{2}} + 1 \right)\frac{4 \pi}{(c h)^3} \int_{\epsilon_0}^{\epsilon_{\mathrm{F}}}d\epsilon \epsilon \sqrt{ \epsilon^2 -\epsilon_0^2} \\ &= \frac{8 \pi}{(c h)^3}\frac{1}{{3}} {\left.{{\left( x^2 - \epsilon_0^2 \right)^{3/2}}}\right\vert}_{{\epsilon_0}}^{{\epsilon_{\mathrm{F}}}} \\ &= \frac{8 \pi}{3 (c h)^3}\left( \epsilon_{\mathrm{F}}^2 - \epsilon_0^2 \right)^{3/2}.\end{aligned} \hspace{\stretch{1}}(1.0.2)

Solving for $\epsilon_{\mathrm{F}}/\epsilon_0$ we have

\begin{aligned}\frac{\epsilon_{\mathrm{F}}}{\epsilon_0} =\sqrt{\left( \frac{3 (c h)^3 n}{8 \pi \epsilon_0^3} \right)^{2/3}+ 1}.\end{aligned} \hspace{\stretch{1}}(1.0.2)

We’ll see the constant factor above a number of times below and designate it

\begin{aligned}n_0 = \frac{8 \pi}{3} \left( \frac{\epsilon_0}{c h} \right)^3,\end{aligned} \hspace{\stretch{1}}(1.0.2)

so that the Fermi energy is

\begin{aligned}\frac{\epsilon_{\mathrm{F}}}{\epsilon_0} =\sqrt{\left( \frac{n}{n_0} \right)^{2/3}+ 1}.\end{aligned} \hspace{\stretch{1}}(1.0.2)

For the pressure calculation, let’s suppose that we have a configuration with a plane in the $x,y$ orientation as in fig. 1.1.

Fig 1.1: Pressure against x,y oriented plane

It’s argued in [4] section 6.4 that the pressure for such a configuration is

\begin{aligned}P = n \int p_z u_z f(\mathbf{u}) d^3 \mathbf{u},\end{aligned} \hspace{\stretch{1}}(1.7)

where $n$ is the number density and $f(\mathbf{u})$ is a normalized distribution function for the velocities. The velocity and momentum components are related by the Hamiltonian equations. From the Hamiltonian eq. 1.1 we find \footnote{ Observe that by squaring and summing one can show that this is equivalent to the standard relativisitic momentum $p_x = \frac{m v_x}{\sqrt{ 1 - \mathbf{u}^2/c^2}}$.} (for the x-component which is representative)

\begin{aligned}u_x \\ &= \frac{\partial {\epsilon}}{\partial {p_x}} \\ &= \frac{\partial {}}{\partial {p_x}}\sqrt{(p c)^2 +\epsilon_0^2} \\ &= \frac{ p_x c^2 }{\sqrt{(p c)^2 +\epsilon_0^2}}.\end{aligned} \hspace{\stretch{1}}(1.8)

For $\alpha \in \{1, 2, 3\}$ we can summarize these velocity-momentum relationships as

\begin{aligned}\frac{u_\alpha}{c} = \frac{ c p_\alpha }{ \epsilon }.\end{aligned} \hspace{\stretch{1}}(1.9)

Should we attempt to calculate the pressure with this parameterization of the velocity space we end up with convergence problems, and can’t express the results in terms of $f^+_\nu(z)$. Let’s try instead with a distribution over momentum space

\begin{aligned}P=n \int \frac{(c p_z)^2}{\epsilon} f(c \mathbf{p}) d^3 (c \mathbf{p}).\end{aligned} \hspace{\stretch{1}}(1.10)

Here the momenta have been scaled to have units of energy since we want to express this integral in terms of energy in the end. Our normalized distribution function is

\begin{aligned}f(c \mathbf{p})\propto \frac{\frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }}}{\int \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }} d^3 (c \mathbf{p})},\end{aligned} \hspace{\stretch{1}}(1.11)

but before evaluating anything, we first want to change our integration variable from momentum to energy. In spherical coordinates our volume element takes the form

\begin{aligned}d^3 (c \mathbf{p}) &= 2 \pi (c p)^2 d (c p) \sin\theta d\theta \\ &= 2 \pi (c p)^2 \frac{d (c p)}{d \epsilon} d \epsilon \sin\theta d\theta.\end{aligned} \hspace{\stretch{1}}(1.12)

Implicit derivatives of

\begin{aligned}c^2 p^2 = \epsilon^2 - \epsilon_0^2,\end{aligned} \hspace{\stretch{1}}(1.13)

gives us

\begin{aligned}\frac{d (c p)}{d\epsilon}= \frac{\epsilon}{c p}=\frac{\epsilon}{\sqrt{\epsilon^2 -\epsilon_0^2}}.\end{aligned} \hspace{\stretch{1}}(1.0.14)

Our momentum volume element becomes

\begin{aligned}d^3 (c \mathbf{p}) \\ &= 2 \pi (c p)^2 \frac{\epsilon}{\sqrt{\epsilon^2 - \epsilon_0^2 }}d \epsilon \sin\theta d\theta \\ &= 2 \pi \left( \epsilon^2 - \epsilon_0^2 \right)\frac{\epsilon}{\sqrt{\epsilon^2 - \epsilon_0^2 }}d \epsilon \sin\theta d\theta \\ &= 2 \pi \epsilon \sqrt{ \epsilon^2 - \epsilon_0^2} d \epsilon \sin\theta d\theta.\end{aligned} \hspace{\stretch{1}}(1.0.14)

For our distribution function, we can now write

\begin{aligned}f(c \mathbf{p}) d^3 (c \mathbf{p})= C \frac{\epsilon \sqrt{ \epsilon^2 - \epsilon_0^2} d \epsilon }{ z^{-1} e^{\beta \epsilon} + 1 }\frac{ 2 \pi \sin\theta d\theta }{ 4 \pi \epsilon_0^3 },\end{aligned} \hspace{\stretch{1}}(1.0.14)

where $C$ is determined by the requirement $\int f(c \mathbf{p}) d^3 (c \mathbf{p}) = 1$

\begin{aligned}C^{-1} = \int_{0}^\infty \frac{(y + 1)\sqrt{ (y + 1)^2 - 1} dy }{ z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 }.\end{aligned} \hspace{\stretch{1}}(1.0.14)

The z component of our momentum can be written in spherical coordinates as

\begin{aligned}(c p_z)^2= (c p)^2 \cos^2\theta= \left( \epsilon^2 - \epsilon_0^2 \right)\cos^2\theta,\end{aligned} \hspace{\stretch{1}}(1.0.18)

Noting that

\begin{aligned}\int_0^\pi \cos^2\theta \sin\theta d\theta =-\int_0^\pi \cos^2\theta d(\cos\theta)= \frac{2}{3},\end{aligned} \hspace{\stretch{1}}(1.0.19)

all the bits come together as

\begin{aligned}P &= \frac{C n}{3 \epsilon_0^3 } \int_{\epsilon_0}^\infty\left( \epsilon^2 - \epsilon_0^2 \right)^{3/2} \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1 }} d \epsilon \\ &= \frac{n \epsilon_0}{3} \int_{1}^\infty\left( x^2 - 1 \right)^{3/2} \frac{1}{{ z^{-1} e^{\beta \epsilon_0 x} + 1 }} dx.\end{aligned} \hspace{\stretch{1}}(1.0.19)

Letting $y = x - 1$, this is

\begin{aligned}P= \frac{C n \epsilon_0}{3} \int_{0}^\infty \frac{ \left( (y + 1)^2 - 1 \right)^{3/2} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 } dy.\end{aligned} \hspace{\stretch{1}}(1.0.19)

We could conceivable expand the numerators of each of these integrals in power series, which could then be evaluated as a sum of $f^+_\nu(z e^{-\beta \epsilon_0})$ terms.

Note that above the Fermi energy $n$ also has an integral representation

\begin{aligned}n &= \left(2\left( \frac{1}{{2}} \right) + 1\right)\int_{\epsilon_0}^\infty d\epsilon \mathcal{D}(\epsilon) \frac{1}{{ z^{-1} e^{\beta \epsilon} + 1}} \\ &= \frac{8 \pi}{(c h)^3} \int_{\epsilon_0}^\infty d\epsilon\frac{\epsilon \sqrt{\epsilon^2 - \epsilon_0^2} }{ z^{-1} e^{\beta \epsilon} + 1} \\ &= \frac{8 \pi \epsilon_0^3}{(c h)^3} \int_{0}^\infty dy\frac{(y + 1)\sqrt{(y + 1)^2 - 1} }{ z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1},\end{aligned} \hspace{\stretch{1}}(1.0.19)

or

\begin{aligned}\boxed{n = \frac{3 n_0}{C}.}\end{aligned} \hspace{\stretch{1}}(1.0.23)

Observe that we can use this result to remove the dependence of pressure on this constant $C$

\begin{aligned}\boxed{\frac{P}{n_0 \epsilon_0}= \int_{0}^\infty dy \frac{ \left( (y + 1)^2 - 1 \right)^{3/2} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 }.}\end{aligned} \hspace{\stretch{1}}(1.0.24)

Now for the average energy difference from the rest energy $\epsilon_0$

\begin{aligned}U &= \left\langle{{\epsilon - \epsilon_0}}\right\rangle \\ &= \int_{\epsilon_0}^\infty d\epsilon \mathcal{D}(\epsilon) f(\epsilon) (\epsilon - \epsilon_0) \\ &= \frac{8 \pi V}{(c h)^3}\int_{\epsilon_0}^\infty d\epsilon \frac{ \epsilon(\epsilon - \epsilon_0) \sqrt{ \epsilon^2 - \epsilon_0 } }{ z^{-1} e^{\beta \epsilon} + 1} \\ &= \frac{8 \pi V \epsilon_0^4}{(c h)^3}\int_{0}^\infty dy\frac{ y ( y - 1 ) \sqrt{ (y + 1)^2 - 1 }}{ z^{-1} e^{\beta \epsilon} + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.24)

So the average energy density difference from the rest energy, relative to the rest energy, is

\begin{aligned}\boxed{\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0} =3 n_0 \int_{0}^\infty dy \frac { y (y + 1)\sqrt{(y + 1)^2 - 1} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 }.}\end{aligned} \hspace{\stretch{1}}(1.0.26)

From eq. 1.0.24 and eq. 1.0.26 we have

\begin{aligned}\frac{1}{{n_0}} &= 3 \frac{V \epsilon_0} {\left\langle{{\epsilon - \epsilon_0}}\right\rangle} \int_{0}^\infty \frac { y (y + 1)\sqrt{(y + 1)^2 - 1} dy } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 } \\ &= \frac{\epsilon_0}{P} \int_{0}^\infty \frac{ \left( (y + 1)^2 - 1 \right)^{3/2} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 } dy,\end{aligned} \hspace{\stretch{1}}(1.0.26)

or

\begin{aligned}P V =\frac{U}{3}\frac{ \int_{0}^\infty \frac{ \left( (y + 1)^2 - 1 \right)^{3/2} } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 } dy}{ \int_{0}^\infty \frac { y (y + 1)\sqrt{(y + 1)^2 - 1} dy } { z^{-1} e^{\beta \epsilon_0 (y + 1)} + 1 }}.\end{aligned} \hspace{\stretch{1}}(1.0.26)

This ratio of integrals is supposed to resolve to 1 and 2 in the low and high density limits. To consider this let’s perform one final non-dimensionalization, writing

\begin{aligned}\begin{aligned} \\ x &= \beta \epsilon_0 y \\ \theta &= \frac{1}{{\beta \epsilon_0}} = \frac{k_{\mathrm{B}} T}{\epsilon_0} \\ \bar{\mu} &= \mu - \epsilon_0 \\ \bar{z} &= e^{\beta \bar{\mu}}.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.29)

The density, pressure, and energy take the form

\begin{aligned}\frac{n}{n_0}= 3 \theta\int_{0}^\infty dx\frac{(\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} }{ \bar{z}^{-1} e^{x} + 1}\end{aligned} \hspace{\stretch{1}}(1.0.30a)

\begin{aligned}\frac{P}{n_0 \epsilon_0}= \theta \int_{0}^\infty dx \frac{ \left( (\theta x + 1)^2 - 1 \right)^{3/2} } { \bar{z}^{-1} e^{x} + 1 }\end{aligned} \hspace{\stretch{1}}(1.0.30b)

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} =3 \theta^2 \int_{0}^\infty dx \frac { x (\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} } { \bar{z}^{-1} e^{x} + 1 }.\end{aligned} \hspace{\stretch{1}}(1.0.30c)

We can rewrite the square roots in the number density and energy density expressions by expanding out the completion of the square

\begin{aligned}(1 + \theta x) \sqrt{ (1 + \theta x)^2 - 1}=(1 + \theta x) \sqrt{ 1 + \theta x + 1 }\sqrt{ 1 + \theta x - 1 }= \sqrt{2 \theta} x^{1/2} (1 + \theta x) \sqrt{ 1 + \frac{\theta x}{2}},\end{aligned} \hspace{\stretch{1}}(1.0.30c)

Expanding the distribution about $\bar{z} e^{-x} = 0$, we have

\begin{aligned}\frac{1}{ \bar{z}^{-1} e^{x} + 1}=\frac{\bar{z} e^{-x}}{ 1 + \bar{z} e^{-x}}=z e^{-x} \sum_{s = 0}^\infty (-1)^s \left( \bar{z} e^{-x} \right)^s,\end{aligned} \hspace{\stretch{1}}(1.0.32)

allowing us to write, in the low density limit with respect to $\bar{z}$

\begin{aligned}\frac{n}{n_0}= 3 \sqrt{2}\theta^{3/2} \sum_{s=0}^\infty(-1)^s\bar{z}^{s + 1}\int_{0}^\infty dx x^{1/2}(1 + \theta x) \sqrt{ 1 + \frac{\theta x}{2}} e^{-x(1 + s)} \end{aligned} \hspace{\stretch{1}}(1.0.33a)

\begin{aligned}\frac{P}{n_0 \epsilon_0}= \theta\sum_{s=0}^\infty(-1)^s\bar{z}^{s + 1} \int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2} e^{-x(1 + s)} \end{aligned} \hspace{\stretch{1}}(1.0.33b)

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} =3 \sqrt{2} \theta^{5/2} \sum_{s=0}^\infty(-1)^s\bar{z}^{s + 1} \int_{0}^\infty dx x^{3/2} (1 + \theta x) \sqrt{ 1 + \frac{\theta x}{2}} e^{-x(1 + s)} .\end{aligned} \hspace{\stretch{1}}(1.0.33c)

Low density result

An exact integration of the various integrals above is possible in terms of special functions. However, that attempt (included below) introduced an erroneous extra factor of $\theta$. Given that this end result was obtained by tossing all but the lowest order terms in $\theta$ and $\bar{z}$, let’s try that right from the get go.

For the pressure we have an integrand containing a factor

\begin{aligned}\left( (\theta x + 1)^2 -1 \right)^{3/2}&= \left( \theta x + 1 - 1 \right)^{3/2}\left( \theta x + 1 + 1 \right)^{3/2} \\ &= \theta^{3/2} x^{3/2} 2^{3/2} \left( 1 + \frac{\theta x}{2} \right)^{3/2} \\ &= 2 \sqrt{2} \theta^{3/2} x^{3/2} \left( 1 + \frac{\theta x}{2} \right)^{3/2}\approx2 \sqrt{2} \theta^{3/2} x^{3/2} \end{aligned} \hspace{\stretch{1}}(1.0.33c)

Our pressure, to lowest order in $\theta$ and $\bar{z}$ is then

\begin{aligned}\frac{P}{\epsilon_0 n_0} = 2 \sqrt{2} \theta^{5/2} \bar{z} \int_0^\infty x^{3/2} e^{-x} dx= 2 \sqrt{2} \theta^{5/2} \bar{z} \Gamma(5/2).\end{aligned} \hspace{\stretch{1}}(1.0.33c)

Our energy density to lowest order in $\theta$ and $\bar{z}$ from eq. 1.0.33c is

\begin{aligned}\frac{U}{V \epsilon_0 n_0} &= 3 \sqrt{2} \theta^{5/2} \bar{z} \int_{0}^\infty dx x^{3/2} e^{-x} \\ &= 3 \sqrt{2} \theta^{5/2} \bar{z} \Gamma(5/2).\end{aligned} \hspace{\stretch{1}}(1.0.33c)

Comparing these, we have

\begin{aligned}\frac{1}{{\epsilon_0 n_0\sqrt{2} \theta^{5/2} \bar{z} \Gamma(5/2)}} &= 3 \frac{V}{U} \\ &= \frac{2}{P},\end{aligned} \hspace{\stretch{1}}(1.0.37)

or in this low density limit

\begin{aligned}\boxed{P V = \frac{2}{3} U.}\end{aligned} \hspace{\stretch{1}}(1.0.38)

High density limit

For the high density limit write $\bar{z} = e^y$, so that the distribution takes the form

\begin{aligned}f(\bar{z}) &= \frac{1}{ \bar{z}^{-1} e^{x} + 1} \\ &= \frac{1}{ e^{x - y} + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.39)

This can be approximated by a step function, so that

\begin{aligned}\frac{P}{n_0 \epsilon_0}\approx \int_{0}^y \theta dx\left( (\theta x + 1)^2 - 1 \right)^{3/2} \end{aligned} \hspace{\stretch{1}}(1.0.40a)

\begin{aligned}\frac{U}{V \epsilon_0 n_0} \approx3 \int_{0}^\infty \theta dx \theta x (\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} \end{aligned} \hspace{\stretch{1}}(1.0.40b)

With a change of variables $u = \theta x + 1$, we have

\begin{aligned}\begin{aligned}\frac{P}{n_0 \epsilon_0} &\approx \int_{1}^{\theta y + 1x} du\left( u^2 - 1 \right)^{3/2} \\ &=\frac{1}{8} \left((2 \theta y (\theta y+2)-3) \sqrt{\theta y (\theta y+2)} (\theta y+1)+3 \ln \left(\theta y+\sqrt{\theta y (\theta y+2)}+1\right)\right) \\ &\approx\frac{1}{4} \left( \theta \ln \bar{z} \right)^4\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.41a)

\begin{aligned}\begin{aligned}\frac{U}{V \epsilon_0 n_0} &\approx3 \int_{1}^{\theta y + 1x} (u^2 - u)\sqrt{u^2 - 1} \\ &=\frac{3}{24} \left(\sqrt{\theta y (\theta y+2)} (\theta y (2 \theta y (3 \theta y+5)-1)+3)-3 \left(\ln \left(\theta y+\sqrt{\theta y (\theta y+2)}+1\right)\right)\right) \\ &\approx\frac{3}{4} \left( \theta \ln \bar{z} \right)^4\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.41b)

Comparing both we have

\begin{aligned}\frac{4}{\epsilon_0 n_0 \left( \theta \ln \bar{z} \right) } = \frac{1}{{P}} = \frac{3 V}{U},\end{aligned} \hspace{\stretch{1}}(1.0.42)

or

\begin{aligned}\boxed{P V = \frac{1}{{3}} U.}\end{aligned} \hspace{\stretch{1}}(1.0.43)

\begin{aligned}{\left.{{\epsilon_{\mathrm{F}}}}\right\vert}_{{n = 1/(0.01)^3}} = 6.12402 \times 10^{-35} \text{J} \times 6.24150934 \times 10^{18} \frac{\text{eV}}{\text{J}} = 3.82231 \times 10^{-16} \text{eV}\end{aligned} \hspace{\stretch{1}}(1.0.43)

Wow. That’s pretty low!

Pressure integral

Of these the pressure integral is yields directly to Mathematica

\begin{aligned}\begin{aligned} \int_{0}^\infty & dx\left( (\theta x + 1)^2 - 1 \right)^{3/2} e^{-x(1 + s)} \\ &=\frac{3 \theta e^{(s+1)/\theta}}{(s + 1)^2} K_2\left( \frac{s+1}{\theta } \right) \\ &=\frac{3 \sqrt{\frac{\pi }{2}} \theta ^{3/2}}{(s+1)^{5/2}}+\frac{45 \sqrt{\frac{\pi }{2}} \theta ^{5/2}}{8 (s+1)^{7/2}}+\frac{315 \sqrt{\frac{\pi }{2}} \theta ^{7/2}}{128 (s+1)^{9/2}}-\frac{945 \sqrt{\frac{\pi }{2}} \theta ^{9/2}}{1024 (s+1)^{11/2}}+\frac{31185 \sqrt{\frac{\pi }{2}} \theta ^{11/2}}{32768 (s+1)^{13/2}} + \cdots\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.45)

where $K_2(z)$ is a modified Bessel function [5] of the second kind as plotted in fig. 1.2.

Fig 1.2: Modified Bessel function of the second kind

Plugging this into the series for the pressure, we have

\begin{aligned}\frac{P}{n_0 \epsilon_0}= 3 \left( \frac{k_{\mathrm{B}} T}{\epsilon_0} \right)^2\sum_{s=0}^\infty(-1)^s\frac{\left( \bar{z} e^{\epsilon_0/k_{\mathrm{B}} T} \right)^{s + 1}}{(s + 1)^2}K_2\left( (s+1) \epsilon_0/k_{\mathrm{B}} T \right).\end{aligned} \hspace{\stretch{1}}(1.0.46)

Plotting the summands $3 (-1)^s \frac{\theta^2}{(s + 1)^2} \left( \bar{z} e^{ 1/\theta} \right)^{s + 1} K_2\left((s+1)/\theta\right)$ for $\bar{z} = 1$ in fig. 1.4 shows that this mix of exponential Bessel and quadratic terms decreases with $s$.

Plotting this sum in fig. 1.3 numerically to 10 terms, shows that we have a function that appears roughly polynomial in $\bar{z}$ and $\theta$.

Fig 1.3: Pressure to ten terms in z and theta

Fig 1.4: Pressure summands

For small $\bar{z}$ it can be seen graphically that there is very little contribution from anything but the $s = 0$ term of this sum. An expansion in series for a few terms in $\bar{z}$ and $\theta$ gives us

\begin{aligned}\begin{aligned}\frac{P}{\epsilon_0 n_0}&=\sqrt{\pi} \theta^{5/2} \left(\frac{3 \bar{z}}{\sqrt{2}}-\frac{3 \bar{z}^2}{8}+\frac{\bar{z}^3}{3 \sqrt{6}}-\frac{3 \bar{z}^4}{32 \sqrt{2}}+\frac{3 \bar{z}^5}{25 \sqrt{10}}\right) \\ &+\sqrt{\pi} \theta^{7/2} \left(\frac{45 \bar{z}}{8 \sqrt{2}}\right) -\frac{45 \bar{z}^2}{128}+\frac{5 \bar{z}^3}{24 \sqrt{6}}-\frac{45 \bar{z}^4}{1024 \sqrt{2}}+\frac{9 \bar{z}^5}{200 \sqrt{10}}\\ &+\sqrt{\pi} \theta^{9/2} \left(\frac{315 \bar{z}}{128 \sqrt{2}}-\frac{315 \bar{z}^2}{4096}+\frac{35 \bar{z}^3}{1152 \sqrt{6}}-\frac{315 \bar{z}^4}{65536 \sqrt{2}}+\frac{63 \bar{z}^5}{16000 \sqrt{10}}\right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.47)

This allows a $k_{\mathrm{B}} T \ll m c^2$ and $\bar{z} \ll 1$ approximation of the pressure

\begin{aligned}\frac{P}{\epsilon_0 n_0} = \frac{3}{2} \sqrt{2 \pi} \bar{z} \theta^{5/2}.\end{aligned} \hspace{\stretch{1}}(1.0.48)

Number density integral

For the number density, it appears that we can evaluate the integral using integration from parts applied to eq. 1.0.30.30

\begin{aligned}\frac{n}{n_0}= \theta\int_{0}^\infty dx\frac{3 (\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} }{ \bar{z}^{-1} e^{x} + 1}=\theta\int_{0}^\infty dx\left( \frac{d}{dx} \left( (\theta x + 1)^2 - 1 \right) ^{3/2} \right)\frac{1}{ \bar{z}^{-1} e^{x} + 1}={\left.{{\theta\left( (\theta x + 1)^2 - 1 \right)^{3/2}\frac{1}{ \bar{z}^{-1} e^{x} + 1}}}\right\vert}_{{0}}^{{\infty}}-\theta\int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2}\frac{ -\bar{z}^{-1} e^{x} }{ \left( \bar{z}^{-1} e^{x} + 1 \right)^2}=\theta\int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2}\frac{ \bar{z} e^{-x} }{ \left( 1 + \bar{z} e^{-x} \right)^2}.\end{aligned} \hspace{\stretch{1}}(1.0.48)

Expanding in series, gives us

\begin{aligned}\frac{n}{n_0}=\theta\sum_{s = 0}^\infty\binom{-2}{s}\bar{z}^{s + 1} \int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2} e^{-x(s + 1)}=3 \theta^2\sum_{s = 0}^\infty\binom{-2}{s}\frac{\left( \bar{z} e^{1/\theta} \right)^{s + 1}}{(s + 1)^2}K_2\left( \frac{s+1}{\theta } \right).\end{aligned} \hspace{\stretch{1}}(1.0.48)

Here the binomial coefficient has the meaning given in the definitions of \statmechchapcite{nonIntegralBinomialSeries}, where for negative integral values of $b$ we have

\begin{aligned}\binom{b}{s}\equiv(-1)^s \frac{-b}{-b + s} \binom{-b+s}{-b}.\end{aligned} \hspace{\stretch{1}}(1.0.51)

Expanding in series to a couple of orders in $\theta$ and $\bar{z}$ we have

\begin{aligned}\frac{n}{n_0} = \frac{\sqrt{2 \pi}}{36} \theta^{1/2} \left(\left(2 \sqrt{3} \bar{z} - 9/\sqrt{2} \right) \bar{z} +18 \right) \bar{z}+\frac{5 \sqrt{ 2 \pi}}{576} \theta^{3/2} \left(\left(4 \sqrt{3} \bar{z} - 27/\sqrt{2}\right) \bar{z} +108 \right) \bar{z}+ \cdots\end{aligned} \hspace{\stretch{1}}(1.0.52)

To first order in $\theta$ and $\bar{z}$ this is

\begin{aligned}\frac{n}{n_0} = \frac{1}{{2}} \sqrt{ 2 \pi } \bar{z} \theta^{1/2},\end{aligned} \hspace{\stretch{1}}(1.0.53)

which allows a relation to pressure

\begin{aligned}P V = 3 N (k_{\mathrm{B}} T)^2 /\epsilon_0.\end{aligned} \hspace{\stretch{1}}(1.0.54)

It’s kind of odd seeming that this is quadratic in temperature. Is there an error?

Energy integral

Starting from eq. 1.0.30c and integrating by parts we have

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} &= 3 \theta^2 \int_{0}^\infty dx \frac { x (\theta x + 1)\sqrt{(\theta x + 1)^2 - 1} } { \bar{z}^{-1} e^{x} + 1 } \\ &= -\theta^2 \int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2}\frac{d}{dx} \left( \frac{x} { \bar{z}^{-1} e^{x} + 1 } \right) \\ &= -\theta^2 \int_{0}^\infty dx\left( (\theta x + 1)^2 - 1 \right)^{3/2}\left( \frac{1} { \bar{z}^{-1} e^{x} + 1 } - \frac{x \bar{z}^{-1} e^{x} } { \left( \bar{z}^{-1} e^{x} + 1 \right)^2 } \right) \\ &= \theta^2 \int_{0}^\infty dx \left( (\theta x + 1)^2 - 1 \right)^{3/2} \frac{ (x - 1)\bar{z}^{-1} e^{x} - 1} { \left( \bar{z}^{-1} e^{x} + 1 \right)^2 } \\ &= \theta^2 \int_{0}^\infty dx \left( (\theta x + 1)^2 - 1 \right)^{3/2} \frac{ (x - 1)\bar{z} e^{-x} - \bar{z}^2 e^{-2 x}} { \left( 1 + \bar{z} e^{-x} \right)^2 } \\ &= \theta^2\sum_{s=0}^\infty \binom{-2}{s} \int_{0}^\infty dx \left( (\theta x + 1)^2 - 1 \right)^{3/2} \left( (x - 1)\bar{z} e^{-x} - \bar{z}^2 e^{-2 x} \right) (\bar{z} e^{-x})^s \\ &= \theta^2\sum_{s=0}^\infty \binom{-2}{s} \bar{z}^{s + 1} \int_{0}^\infty dx \left( (\theta x + 1)^2 - 1 \right)^{3/2} \left( (x - 1) e^{-x(s + 1)} - \bar{z} e^{-x(s + 2)} \right).\end{aligned} \hspace{\stretch{1}}(1.0.54)

The integral with the factor of $x$ doesn’t have a nice closed form as before (if you consider the $K_2$ a nice closed form), but instead evaluates to a confluent hypergeometric function [6]. That integral is

\begin{aligned}\int_0^{\infty } x \left((\theta x+1)^2-1\right)^{3/2} e^{-x (1+s)} dx = \frac{15 \sqrt{\pi } \theta^3 U\left(-\frac{3}{2},-4,\frac{2 (s+1)}{\theta }\right)}{8 (s+1)^5},\end{aligned} \hspace{\stretch{1}}(1.0.54)

and looks like fig. 1.5. Series expansion shows that this hypergeometricU function has a $\theta^{3/2}$ singularity at the origin

Fig 1.5: Plot of HypergeometricU, and with theta^5 scaling

\begin{aligned}U\left(-\frac{3}{2},-4,\frac{2 (s+1)}{\theta }\right)=\frac{2 \sqrt{2} \sqrt{s+1} s+2 \sqrt{2} \sqrt{s+1}}{\theta^{3/2}}+\frac{21 \sqrt{s+1}}{2 \sqrt{2} \sqrt{\theta }}+ \cdots\end{aligned} \hspace{\stretch{1}}(1.57)

so our multiplication by $\theta^5$ brings us to zero as seen in the plot. Evaluating the complete integral yields the unholy mess

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} &= \sum_{s=0}^\infty \theta^2 (-1)^s (s+1) \bar{z}^{s+1} \Bigl( \\ &\frac{105 \sqrt{\pi } \theta^3 U\left(-\frac{1}{2},-4,\frac{2 (s+1)}{\theta }\right)}{16 (s+1)^5} \\ &- \frac{3 \sqrt{\pi } \theta^2 U\left(-\frac{1}{2},-2,\frac{2 (s+1)}{\theta }\right)}{2 (s+1)^3} \\ &- \frac{3 \sqrt{\pi } \theta^2 \bar{z} U\left(-\frac{1}{2},-2,\frac{2 (s+2)}{\theta }\right)}{2 (s+2)^3} \\ &+\frac{(\theta -2) (-3 \theta +2 s+2) e^{\frac{s+1}{\theta }} K_2\left(\frac{s+1}{\theta }\right)}{\theta (s+1)^2} \\ &-\frac{2 (\theta -2) e^{\frac{s+1}{\theta }} K_1\left(\frac{s+1}{\theta }\right)}{\theta (s+1)} \\ &+\frac{\bar{z} (-3 \theta +2 s+4) e^{\frac{s+2}{\theta }} K_2\left(\frac{s+2}{\theta }\right)}{(s+2)^2} \\ &-\frac{2 \bar{z} e^{\frac{s+2}{\theta }} K_1\left(\frac{s+2}{\theta }\right)}{s+2} \Bigr),\end{aligned} \hspace{\stretch{1}}(1.58)

to first order in $\bar{z}$ and $\theta$ this is

\begin{aligned}\frac{\left\langle{{\epsilon - \epsilon_0}}\right\rangle}{V \epsilon_0 n_0} =\frac{9}{4} \sqrt{2 \pi} \bar{z} \theta^{7/2}.\end{aligned} \hspace{\stretch{1}}(1.59)

Comparing pressure and energy we have for low densities (where $\bar{z} \approx 0$)

\begin{aligned}\frac{1}{{\epsilon_0 n_0 \sqrt{2 \pi} \bar{z} \theta^{5/2}}} = \frac{3}{2} \frac{1}{{P}} = \frac{9}{4} \theta \frac{V}{U},\end{aligned} \hspace{\stretch{1}}(1.0.60)

or

\begin{aligned}\theta P V = \frac{2}{3} U.\end{aligned} \hspace{\stretch{1}}(1.0.61)

It appears that I’ve picked up an extra factor of $\theta$ somewhere, but at least I’ve got the $2/3$ low density expression. Given that I’ve Taylor expanded everything anyways around $\bar{z}$ and $\theta$ this could likely have been done right from the get go, instead of dragging along the messy geometric integrals. Reworking this part of this problem like that was done above.

References

[1] Kerson Huang. Introduction to statistical physics. CRC Press, 2001.

[2] Peeter Joot. Basic statistical mechanics., chapter {Non integral binomial coefficient}. \natexlab{a}. URL http://sites.google.com/site/peeterjoot2/math2013/phy452.pdf.

[3] Peeter Joot. Basic statistical mechanics., chapter {Relativisitic density of states}. \natexlab{b}. URL http://sites.google.com/site/peeterjoot2/math2013/phy452.pdf.

[4] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996.

[5] Wolfram. BesselK, \natexlab{a}. URL http://reference.wolfram.com/mathematica/ref/BesselK.html. [Online; accessed 11-April-2013].

[6] Wolfram. HyperGeometricU, \natexlab{b}. URL http://reference.wolfram.com/mathematica/ref/HypergeometricU.html. [Online; accessed 17-April-2013].

An updated compilation of notes, for ‘PHY452H1S Basic Statistical Mechanics’, Taught by Prof. Arun Paramekanti

Posted by peeterjoot on March 27, 2013

Here’s my second update of my notes compilation for this course, including all of the following:

March 27, 2013 Fermi gas

March 26, 2013 Fermi gas thermodynamics

March 26, 2013 Fermi gas thermodynamics

March 23, 2013 Relativisitic generalization of statistical mechanics

March 21, 2013 Kittel Zipper problem

March 18, 2013 Pathria chapter 4 diatomic molecule problem

March 17, 2013 Gibbs sum for a two level system

March 16, 2013 open system variance of N

March 16, 2013 probability forms of entropy

March 14, 2013 Grand Canonical/Fermion-Bosons

March 13, 2013 Quantum anharmonic oscillator

March 12, 2013 Grand canonical ensemble

March 11, 2013 Heat capacity of perturbed harmonic oscillator

March 10, 2013 Langevin small approximation

March 10, 2013 Addition of two one half spins

March 10, 2013 Midterm II reflection

March 07, 2013 Thermodynamic identities

March 06, 2013 Temperature

March 05, 2013 Interacting spin

plus everything detailed in the description of my first update and before.

PHY452H1S Basic Statistical Mechanics. Lecture 16: Fermi gas. Taught by Prof. Arun Paramekanti

Posted by peeterjoot on March 27, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Disclaimer

Peeter’s lecture notes from class. May not be entirely coherent.

Fermi gas

Review

Continuing a discussion of [1] section 8.1 content.

We found

\begin{aligned}n_{\mathbf{k}} = \frac{1}{{e^{\beta(\epsilon_k - \mu)} + 1}}\end{aligned} \hspace{\stretch{1}}(1.2.1)

With no spin

\begin{aligned}\int n_\mathbf{k} \times \frac{d^3 k}{(2\pi)^3} = \rho\end{aligned} \hspace{\stretch{1}}(1.2.2)

Fig 1.1: Occupancy at low temperature limit

Fig 1.2: Volume integral over momentum up to Fermi energy limit

\begin{aligned}\epsilon_{\mathrm{F}} = \frac{\hbar^2 k_{\mathrm{F}}^2}{2m}\end{aligned} \hspace{\stretch{1}}(1.2.3)

gives

\begin{aligned}k_{\mathrm{F}} = (6 \pi^2 \rho)^{1/3}\end{aligned} \hspace{\stretch{1}}(1.2.4)

\begin{aligned}\sum_\mathbf{k} n_\mathbf{k} = N\end{aligned} \hspace{\stretch{1}}(1.2.5)

\begin{aligned}\mathbf{k} = \frac{2\pi}{L}(n_x, n_y, n_z)\end{aligned} \hspace{\stretch{1}}(1.2.6)

This is for periodic boundary conditions \footnote{I filled in details in the last lecture using a particle in a box, whereas this periodic condition was intended. We see that both achieve the same result}, where

\begin{aligned}\Psi(x + L) = \Psi(x)\end{aligned} \hspace{\stretch{1}}(1.2.7)

Moving on

\begin{aligned}\sum_{k_x} n(\mathbf{k}) = \sum_{p_x} \Delta p_x n(\mathbf{k})\end{aligned} \hspace{\stretch{1}}(1.2.8)

with

\begin{aligned}\Delta k_x = \frac{2 \pi}{L} \Delta p_x\end{aligned} \hspace{\stretch{1}}(1.2.9)

this gives

\begin{aligned}\sum_{k_x} n(\mathbf{k}) = \sum_{n_x} \frac{L}{2\pi} \Delta k_x \rightarrow \frac{L}{2\pi} \int d k_x\end{aligned} \hspace{\stretch{1}}(1.2.10)

Over all dimensions

\begin{aligned}\sum_{\mathbf{k}} n_\mathbf{k} = \left( \frac{L}{2\pi} \right)^3 \left( \int d^3 \mathbf{k} \right)n(\mathbf{k})=N\end{aligned} \hspace{\stretch{1}}(1.2.11)

so that

\begin{aligned}\rho = \int \frac{d^3 \mathbf{k}}{(2 \pi)^3}\end{aligned} \hspace{\stretch{1}}(1.2.12)

Again

\begin{aligned}k_{\mathrm{F}} = (6 \pi^2 \rho)^{1/3}\end{aligned} \hspace{\stretch{1}}(1.2.13)

Example: Spin considerations

{example:basicStatMechLecture16:1}{

\begin{aligned}\sum_{\mathbf{k}, m_s} = N\end{aligned} \hspace{\stretch{1}}(1.2.14)

\begin{aligned}\sum_{\mathbf{k}, m_s} \frac{1}{{e^{\beta(\epsilon_k - \mu)} + 1}} = (2 S + 1)\left( \int \frac{d^3 \mathbf{k}}{(2 \pi)^3} n(\mathbf{k}) \right)L^3\end{aligned} \hspace{\stretch{1}}(1.2.15)

This gives us

\begin{aligned}k_{\mathrm{F}} = \left( \frac{ 6 \pi^2 \rho }{2 S + 1} \right)^{1/3}\end{aligned} \hspace{\stretch{1}}(1.2.16)

and again

\begin{aligned}\epsilon_{\mathrm{F}} = \frac{\hbar^2 k_{\mathrm{F}}^2}{2m}\end{aligned} \hspace{\stretch{1}}(1.2.17)

}

High Temperatures

Now we want to look at the at higher temperature range, where the occupancy may look like fig. 1.3

Fig 1.3: Occupancy at higher temperatures

\begin{aligned}\mu(T = 0) = \epsilon_{\mathrm{F}}\end{aligned} \hspace{\stretch{1}}(1.2.18)

\begin{aligned}\mu(T \rightarrow \infty) \rightarrow - \infty\end{aligned} \hspace{\stretch{1}}(1.2.19)

so that for large $T$ we have

\begin{aligned}\frac{1}{{e^{\beta(\epsilon_k - \mu)} + 1}} \rightarrow e^{-\beta(\epsilon_k - \mu)}\end{aligned} \hspace{\stretch{1}}(1.2.20)

\begin{aligned}\rho &= \int \frac{d^3 \mathbf{k}}{(2 \pi)^3} e^{\beta \mu} e^{-\beta \epsilon_k} \\ &= e^{\beta \mu} \int \frac{d^3 \mathbf{k}}{(2 \pi)^3} e^{-\beta \epsilon_k} \\ &= e^{\beta \mu} \int dk \frac{4 \pi k^2}{(2 \pi)^3} e^{-\beta \hbar^2 k^2/2m}.\end{aligned} \hspace{\stretch{1}}(1.2.21)

Mathematica (or integration by parts) tells us that

\begin{aligned}\frac{1}{{(2 \pi)^3}} \int 4 \pi^2 k^2 dk e^{-a k^2} = \frac{1}{{(4 \pi a )^{3/2}}},\end{aligned} \hspace{\stretch{1}}(1.2.22)

so we have

\begin{aligned}\rho &= e^{\beta \mu} \left( \frac{2m}{ 4 \pi \beta \hbar^2} \right)^{3/2} \\ &= e^{\beta \mu} \left( \frac{2 m k_{\mathrm{B}} T 4 \pi^2 }{ 4 \pi h^2} \right)^{3/2} \\ &= e^{\beta \mu} \left( \frac{2 m k_{\mathrm{B}} T \pi }{ h^2} \right)^{3/2}\end{aligned} \hspace{\stretch{1}}(1.2.23)

Introducing $\lambda$ for the thermal de Broglie wavelength, $\lambda^3 \sim T^{-3/2}$

\begin{aligned}\lambda \equiv \frac{h}{\sqrt{2 \pi m k_{\mathrm{B}} T}},\end{aligned} \hspace{\stretch{1}}(1.2.24)

we have

\begin{aligned}\rho = e^{\beta \mu} \frac{1}{{\lambda^3}}.\end{aligned} \hspace{\stretch{1}}(1.2.25)

Does it make any sense to have density as a function of temperature? An inappropriately extended to low temperatures plot of the density is found in fig. 1.4 for a few arbitrarily chosen numerical values of the chemical potential $\mu$, where we see that it drops to zero with temperature. I suppose that makes sense if we are not holding volume constant.

Fig 1.4: Density as a function of temperature

We can write

\begin{aligned}\boxed{e^{\beta \mu} = \left( \rho \lambda^3 \right)}\end{aligned} \hspace{\stretch{1}}(1.2.26)

\begin{aligned}\frac{\mu}{k_{\mathrm{B}} T} = \ln \left( \rho \lambda^3 \right)\sim -\frac{3}{2} \ln T\end{aligned} \hspace{\stretch{1}}(1.2.27)

or (taking $\rho$ (and/or volume?) as a constant) we have for large temperatures

\begin{aligned}\mu \propto -T \ln T\end{aligned} \hspace{\stretch{1}}(1.2.28)

The chemical potential is plotted in fig. 1.5, whereas this $- k_{\mathrm{B}} T \ln k_{\mathrm{B}} T$ function is plotted in fig. 1.6. The contributions to $\mu$ from the $k_{\mathrm{B}} T \ln (\rho h^3 (2 \pi m)^{-3/2})$ term are dropped for the high temperature approximation.

Fig 1.5: Chemical potential over degenerate to classical range

Fig 1.6: High temp approximation of chemical potential, extended back to T = 0

Pressure

\begin{aligned}P = - \frac{\partial {E}}{\partial {V}}\end{aligned} \hspace{\stretch{1}}(1.2.29)

For a classical ideal gas as in fig. 1.7 we have

Fig 1.7: Ideal gas pressure vs volume

\begin{aligned}P = \rho k_{\mathrm{B}} T\end{aligned} \hspace{\stretch{1}}(1.2.30)

For a Fermi gas at $T = 0$ we have

\begin{aligned}E &= \sum_\mathbf{k} \epsilon_k n_k \\ &= \sum_\mathbf{k} \epsilon_k \Theta(\mu_0 - \epsilon_k) \\ &= \frac{V}{(2\pi)^3} \int_{\epsilon_k < \mu_0} \frac{\hbar^2 \mathbf{k}^2}{2 m} d^3 \mathbf{k} \\ &= \frac{V}{(2\pi)^3} \int_0^{k_{\mathrm{F}}} \frac{\hbar^2 \mathbf{k}^2}{2 m} d^3 \mathbf{k} \\ &= \frac{V}{(2\pi)^3} \frac{\hbar^2}{2 m} \int_0^{k_{\mathrm{F}}} k^2 4 \pi k^2 d k\propto k_{\mathrm{F}}^5\end{aligned} \hspace{\stretch{1}}(1.2.31)

Specifically,

\begin{aligned}E(T = 0) = V \times \frac{3}{5} \underbrace{\epsilon_{\mathrm{F}}}_{\sim k_{\mathrm{F}}^2}\underbrace{\rho}_{\sim k_{\mathrm{F}}^3}\end{aligned} \hspace{\stretch{1}}(1.2.32)

or

\begin{aligned}\frac{E}{N} = \frac{3}{5} \epsilon_{\mathrm{F}}\end{aligned} \hspace{\stretch{1}}(1.2.33)

\begin{aligned}E = \frac{3}{5} N \frac{\hbar^2}{2 m} \left( 6 \pi^2 \frac{N}{V} \right)^{2/3} = a V^{-2/3},\end{aligned} \hspace{\stretch{1}}(1.2.34)

so that

\begin{aligned}\frac{\partial {E}}{\partial {V}} = -\frac{2}{3} a V^{-5/3}.\end{aligned} \hspace{\stretch{1}}(1.2.35)

\begin{aligned}P &= -\frac{\partial {E}}{\partial {V}} \\ &= \frac{2}{3} \left( a V^{-2/3} \right)V^{-1} \\ &= \frac{2}{3} \frac{E}{V} \\ &= \frac{2}{3} \left( \frac{3}{5} \epsilon_{\mathrm{F}} \rho \right) \\ &= \frac{2}{5} \epsilon_{\mathrm{F}} \rho.\end{aligned} \hspace{\stretch{1}}(1.2.36)

We see that the pressure ends up deviating from the classical result at low temperatures, as sketched in fig. 1.8. This low temperature limit for the pressure $2 \epsilon_{\mathrm{F}} \rho/5$ is called the degeneracy pressure.

Fig 1.8: Fermi degeneracy pressure

References

[1] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996.

Pathria chapter 4 diatomic molecule problem

Posted by peeterjoot on March 18, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Question: Diatomic molecule ([1] pr 4.7)

Consider a classical system of non-interacting, diatomic molecules enclosed in a box of volume $V$ at temperature $T$. The Hamiltonian of a single molecule is given by

\begin{aligned}H(\mathbf{r}_1, \mathbf{r}_2, \mathbf{p}_1, \mathbf{p}_2) = \frac{1}{{2m}} \left( \mathbf{p}_1^2 + \mathbf{p}_2^2 \right)+\frac{1}{{2}} K \left\lvert {\mathbf{r}_1 - \mathbf{r}_2} \right\rvert^2.\end{aligned} \hspace{\stretch{1}}(1.0.1)

Study the thermodynamics of this system, including the dependence of the quantity $\left\langle{{r_{12}^2}}\right\rangle$ on $T$.

Answer

Partition function
First consider the partition function for a single diatomic pair

\begin{aligned}Z_1 &= \frac{1}{{h^6}} \int d^6 \mathbf{p} d^6 \mathbf{r} e^{-\beta \frac{ \mathbf{p}_1^2 + \mathbf{p}_2^2 }{2m}} e^{-\beta K\frac{ \left\lvert {\mathbf{r}_1 - \mathbf{r}_2} \right\rvert^2 }{2}} \\ &= \frac{1}{{h^6}} \left( \frac{2 \pi m}{\beta} \right)^{6/2}\int d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 e^{-\beta K\frac{ \left\lvert {\mathbf{r}_1 - \mathbf{r}_2} \right\rvert^2 }{2}}\end{aligned} \hspace{\stretch{1}}(1.0.2)

Now we can make a change of variables to simplify the exponential. Let’s write

\begin{aligned}\mathbf{u} = \mathbf{r}_1 - \mathbf{r}_2\end{aligned} \hspace{\stretch{1}}(1.0.3a)

\begin{aligned}\mathbf{v} = \mathbf{r}_2,\end{aligned} \hspace{\stretch{1}}(1.0.3b)

or

\begin{aligned}\mathbf{r}_2 = \mathbf{v}\end{aligned} \hspace{\stretch{1}}(1.0.4a)

\begin{aligned}\mathbf{r}_1=\mathbf{u} + \mathbf{v}.\end{aligned} \hspace{\stretch{1}}(1.0.4b)

Our volume element is

\begin{aligned}d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 = d^3 \mathbf{u} d^3 \mathbf{v} \frac{\partial(\mathbf{r}_1, \mathbf{r}_2)}{\partial(\mathbf{u}, \mathbf{v})}.\end{aligned} \hspace{\stretch{1}}(1.0.5)

It wasn’t obvious to me that this change of variables preserves the volume element, but a quick Jacobian calculation shows this to be the case

\begin{aligned}\frac{\partial(\mathbf{r}_1, \mathbf{r}_2)}{\partial(\mathbf{u}, \mathbf{v})} &= \begin{vmatrix}\partial r_{11}/\partial u_1 & \partial r_{11}/\partial u_2 &\partial r_{11}/\partial u_3 &\partial r_{11}/\partial v_1 &\partial r_{11}/\partial v_2 &\partial r_{11}/\partial v_3 \\ \partial r_{12}/\partial u_1 & \partial r_{12}/\partial u_2 &\partial r_{12}/\partial u_3 &\partial r_{12}/\partial v_1 &\partial r_{12}/\partial v_2 &\partial r_{12}/\partial v_3 \\ \partial r_{13}/\partial u_1 & \partial r_{13}/\partial u_2 &\partial r_{13}/\partial u_3 &\partial r_{13}/\partial v_1 &\partial r_{13}/\partial v_2 &\partial r_{13}/\partial v_3 \\ \partial r_{21}/\partial u_1 & \partial r_{21}/\partial u_2 &\partial r_{21}/\partial u_3 &\partial r_{21}/\partial v_1 &\partial r_{21}/\partial v_2 &\partial r_{21}/\partial v_3 \\ \partial r_{22}/\partial u_1 & \partial r_{22}/\partial u_2 &\partial r_{22}/\partial u_3 &\partial r_{22}/\partial v_1 &\partial r_{22}/\partial v_2 &\partial r_{22}/\partial v_3 \\ \partial r_{23}/\partial u_1 & \partial r_{23}/\partial u_2 &\partial r_{23}/\partial u_3 &\partial r_{23}/\partial v_1 &\partial r_{23}/\partial v_2 &\partial r_{23}/\partial v_3 \end{vmatrix} \\ &= \begin{vmatrix}\partial r_{11}/\partial u_1 & \partial r_{11}/\partial u_2 &\partial r_{11}/\partial u_3 &\partial r_{11}/\partial v_1 &\partial r_{11}/\partial v_2 &\partial r_{11}/\partial v_3 \\ \partial r_{12}/\partial u_1 & \partial r_{12}/\partial u_2 &\partial r_{12}/\partial u_3 &\partial r_{12}/\partial v_1 &\partial r_{12}/\partial v_2 &\partial r_{12}/\partial v_3 \\ \partial r_{13}/\partial u_1 & \partial r_{13}/\partial u_2 &\partial r_{13}/\partial u_3 &\partial r_{13}/\partial v_1 &\partial r_{13}/\partial v_2 &\partial r_{13}/\partial v_3 \\ 0 & 0 & 0 &\partial r_{21}/\partial v_1 &\partial r_{21}/\partial v_2 &\partial r_{21}/\partial v_3 \\ 0 & 0 & 0 &\partial r_{22}/\partial v_1 &\partial r_{22}/\partial v_2 &\partial r_{22}/\partial v_3 \\ 0 & 0 & 0 &\partial r_{23}/\partial v_1 &\partial r_{23}/\partial v_2 &\partial r_{23}/\partial v_3 \end{vmatrix} \\ &= 1.\end{aligned} \hspace{\stretch{1}}(1.0.6)

Our remaining integral can now be evaluated

\begin{aligned}\int d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 e^{-\beta K\frac{ \left\lvert {\mathbf{r}_1 - \mathbf{r}_2} \right\rvert^2 }{2}} &= \int d^3 \mathbf{u} d^3 \mathbf{v} e^{-\beta K \left\lvert {\mathbf{u}} \right\rvert^2 /2 } \\ &= V \int d^3 \mathbf{u} e^{-\beta K \left\lvert {\mathbf{u}} \right\rvert^2 /2 } \\ &= V \int d^3 \mathbf{u} e^{-\beta K \left\lvert {\mathbf{u}} \right\rvert^2 /2 } \\ &= V \left( \frac{ 2 \pi }{ K \beta } \right)^{3/2}.\end{aligned} \hspace{\stretch{1}}(1.0.7)

Our partition function is now completely evaluated

\begin{aligned}Z_1 = V\frac{1}{{h^6}} \left( \frac{2 \pi m}{\beta} \right)^{3}\left( \frac{ 2 \pi }{ K \beta } \right)^{3/2}.\end{aligned} \hspace{\stretch{1}}(1.0.8)

As a function of $V$ and $T$ as in the text, we write

\begin{aligned}Z_1 = V f(T)\end{aligned} \hspace{\stretch{1}}(1.0.9a)

\begin{aligned}f(T) = \left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{9/2}.\end{aligned} \hspace{\stretch{1}}(1.0.9b)

Gibbs sum

Our Gibbs sum, summing over the number of molecules (not atoms), is

\begin{aligned}Z_{\mathrm{G}} &= \sum_{N_r = 0}^\infty \frac{z^{N_r}}{N_r!} Z_1^{N_r} \\ &= e^{ z V f(T) },\end{aligned} \hspace{\stretch{1}}(1.0.10)

or

\begin{aligned}q &= \ln Z_{\mathrm{G}} \\ &= z V f(T) \\ &= P V \beta.\end{aligned} \hspace{\stretch{1}}(1.0.11)

The fact that we can sum this as an exponential series so nicely looks like it’s one of the main advantages to this grand partition function (Gibbs sum). We can avoid any of the large $N!$ approximations that we have to use when the number of particles is explicitly fixed.

Pressure

The pressure follows

\begin{aligned}P &= z f(T) k_{\mathrm{B}} T \\ &= e^{\mu/k_{\mathrm{B}} T}\left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{11/2}.\end{aligned} \hspace{\stretch{1}}(1.0.12)

Average energy

\begin{aligned}\left\langle{{H}}\right\rangle &= -\frac{\partial {q}}{\partial {\beta}} \\ &= - z V \frac{9}{2} \frac{f(T)}{T} \frac{\partial {T}}{\partial {\beta}} \\ &= z V \frac{9}{2} \frac{f(T)}{T^3} \frac{1}{{k_{\mathrm{B}}}},\end{aligned} \hspace{\stretch{1}}(1.0.13)

or

\begin{aligned}\left\langle{{H}}\right\rangle = e^{\mu/k_{\mathrm{B}} T} V \frac{9}{2} k_{\mathrm{B}}^2 \left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{3/2}.\end{aligned} \hspace{\stretch{1}}(1.0.14)

Average occupancy

\begin{aligned}\left\langle{{N}}\right\rangle &= z \frac{\partial {}}{\partial {z}} \ln Z_{\mathrm{G}} \\ &= z \frac{\partial {}}{\partial {z}} \left( z V f(T) \right) \\ &= z V f(T)\end{aligned} \hspace{\stretch{1}}(1.0.15)

but this is just $q$, or

\begin{aligned}\left\langle{{N}}\right\rangle &= e^{\mu/k_{\mathrm{B}} T} V\left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{9/2}.\end{aligned} \hspace{\stretch{1}}(1.0.16)

Free energy

\begin{aligned}F &= - k_{\mathrm{B}} T \ln \frac{ Z_{\mathrm{G}} }{z^N} \\ &= - k_{\mathrm{B}} T \left( q - N \ln z \right) \\ &= N k_{\mathrm{B}} T \beta \mu - k_{\mathrm{B}} T q \\ &= z V f(T) \mu - k_{\mathrm{B}} T z V f(T) \\ &= z V f(T) \left( \mu - k_{\mathrm{B}} T \right)\end{aligned} \hspace{\stretch{1}}(1.0.17)

\begin{aligned}F = e^{\mu/k_{\mathrm{B}} T} V \left( \mu - k_{\mathrm{B}} T \right)\left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( k_{\mathrm{B}} T \right)^{9/2}.\end{aligned} \hspace{\stretch{1}}(1.0.18)

Entropy

\begin{aligned}S &= \frac{U - F}{T} \\ &= \frac{V}{T} e^{\mu/k_{\mathrm{B}} T} \left( k_{\mathrm{B}} T \right)^{3/2}\left( \frac{m }{h^2 } \sqrt{\frac{(2\pi)^3}{K}} \right)^3\left( \frac{9}{2} k_{\mathrm{B}}^2 - \left( \mu - k_{\mathrm{B}} T \right) \left( k_{\mathrm{B}} T \right)^3 \right).\end{aligned} \hspace{\stretch{1}}(1.0.19)

Expectation of atomic separation

The momentum portions of the average will just cancel out, leaving just

\begin{aligned}\left\langle{r_{12}^2}\right\rangle &= \frac{\int d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 \left( \mathbf{r}_1 - \mathbf{r}_2 \right)^2 e^{-\beta K \left( \mathbf{r}_1 - \mathbf{r}_2 \right)^2 /2 }}{\int d^3 \mathbf{r}_1 d^3 \mathbf{r}_2 e^{-\beta K \left( \mathbf{r}_1 - \mathbf{r}_2 \right)^2 /2 }} \\ &= \frac{ \int d^3 \mathbf{u} \mathbf{u}^2 e^{-\beta K \mathbf{u}^2 /2 }}{\int d^3 \mathbf{u} e^{-\beta K \mathbf{u}^2 /2 }} \\ &= \frac{\int da db dc \left( a^2 + b^2 + c^2 \right) e^{-\beta K \left( a^2 + b^2 + c^2 \right) /2}}{\int e^{-\beta K \left( a^2 + b^2 + c^2 \right)/2}} \\ &= 3 \frac{\int da a^2 e^{-\beta K a^2/2}\int db dc e^{-\beta K \left( b^2 + c^2 \right) /2}}{\int e^{-\beta K \left( a^2 + b^2 + c^2 \right)/2 }} \\ &= 3 \frac{\int da a^2 e^{-\beta K a^2/2}}{\int e^{-\beta K a^2/2}}\end{aligned} \hspace{\stretch{1}}(1.0.20)

Expanding the numerator by parts we have

\begin{aligned}\int da a^2 e^{-\beta K a^2/2} \\ &= \int a d\frac{ e^{-\beta K a^2/2}}{- 2 \beta K/2} \\ &= \frac{1}{\beta K}\int e^{-\beta K a^2/2}.\end{aligned} \hspace{\stretch{1}}(1.0.21)

This gives us

\begin{aligned}\boxed{\left\langle r_{12}^2 \right\rangle = \frac{3}{\beta K} = \frac{3 k_{\mathrm{B}} T}{K}.}\end{aligned} \hspace{\stretch{1}}(1.0.22)

References

[1] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996.

Thermodynamic identities

Posted by peeterjoot on March 7, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Impressed with the clarity of Baez’s entropic force discussion on differential forms [1], let’s use that methodology to find all the possible identities that we can get from the thermodynamic identity (for now assuming $N$ is fixed, ignoring the chemical potential.)

This isn’t actually that much work to do, since a bit of editor regular expression magic can do most of the work.

Our starting point is the thermodynamic identity

\begin{aligned}dU = d Q + d W = T dS - P dV,\end{aligned} \hspace{\stretch{1}}(1.0.1)

or

\begin{aligned}0 = dU - T dS + P dV.\end{aligned} \hspace{\stretch{1}}(1.0.2)

It’s quite likely that many of the identities that can be obtained will be useful, but this should at least provide a handy reference of possible conversions.

Differentials in $P, V$

This first case illustrates the method.

\begin{aligned}0 &= dU - T dS + P dV \\ &= \left( \frac{\partial {U}}{\partial {P}} \right)_{V} dP +\left( \frac{\partial {U}}{\partial {V}} \right)_{P} dV- T\left( \left( \frac{\partial {S}}{\partial {P}} \right)_{V} dP + \left( \frac{\partial {S}}{\partial {V}} \right)_{P} dV \right)+ P dV \\ &= dP \left( \left( \frac{\partial {U}}{\partial {P}} \right)_{V} - T \left( \frac{\partial {S}}{\partial {P}} \right)_{V} \right)+dV \left( \left( \frac{\partial {U}}{\partial {V}} \right)_{P} - T \left( \frac{\partial {S}}{\partial {V}} \right)_{P} + P \right).\end{aligned} \hspace{\stretch{1}}(1.0.3)

Taking wedge products with $dV$ and $dP$ respectively, we form two two forms

\begin{aligned}0 = dP \wedge dV \left( \left( \frac{\partial {U}}{\partial {P}} \right)_{V} - T \left( \frac{\partial {S}}{\partial {P}} \right)_{V} \right)\end{aligned} \hspace{\stretch{1}}(1.0.4a)

\begin{aligned}0 = dV \wedge dP \left( \left( \frac{\partial {U}}{\partial {V}} \right)_{P} - T \left( \frac{\partial {S}}{\partial {V}} \right)_{P} + P \right).\end{aligned} \hspace{\stretch{1}}(1.0.4b)

Since these must both be zero we find

\begin{aligned}\left( \frac{\partial {U}}{\partial {P}} \right)_{V} = T \left( \frac{\partial {S}}{\partial {P}} \right)_{V}\end{aligned} \hspace{\stretch{1}}(1.0.5a)

\begin{aligned}P =-\left( \frac{\partial {U}}{\partial {V}} \right)_{P}- T \left( \frac{\partial {S}}{\partial {V}} \right)_{P}.\end{aligned} \hspace{\stretch{1}}(1.0.5b)

Differentials in $P, T$

\begin{aligned}0 &= dU - T dS + P dV \\ &= \left( \frac{\partial {U}}{\partial {P}} \right)_{T} dP + \left( \frac{\partial {U}}{\partial {T}} \right)_{P} dT-T \left( \left( \frac{\partial {S}}{\partial {P}} \right)_{T} dP + \left( \frac{\partial {S}}{\partial {T}} \right)_{P} dT \right)+\left( \frac{\partial {V}}{\partial {P}} \right)_{T} dP + \left( \frac{\partial {V}}{\partial {T}} \right)_{P} dT,\end{aligned} \hspace{\stretch{1}}(1.0.6)

or

\begin{aligned}0 = \left( \frac{\partial {U}}{\partial {P}} \right)_{T} -T \left( \frac{\partial {S}}{\partial {P}} \right)_{T} + \left( \frac{\partial {V}}{\partial {P}} \right)_{T}\end{aligned} \hspace{\stretch{1}}(1.0.7a)

\begin{aligned}0 = \left( \frac{\partial {U}}{\partial {T}} \right)_{P} -T \left( \frac{\partial {S}}{\partial {T}} \right)_{P} + \left( \frac{\partial {V}}{\partial {T}} \right)_{P}.\end{aligned} \hspace{\stretch{1}}(1.0.7b)

Differentials in $P, S$

\begin{aligned}0 &= dU - T dS + P dV \\ &= \left( \frac{\partial {U}}{\partial {P}} \right)_{S} dP + \left( \frac{\partial {U}}{\partial {S}} \right)_{P} dS- T dS+ P \left( \left( \frac{\partial {V}}{\partial {P}} \right)_{S} dP + \left( \frac{\partial {V}}{\partial {S}} \right)_{P} dS \right),\end{aligned} \hspace{\stretch{1}}(1.0.8)

or

\begin{aligned}\left( \frac{\partial {U}}{\partial {P}} \right)_{S} = -P \left( \frac{\partial {V}}{\partial {P}} \right)_{S}\end{aligned} \hspace{\stretch{1}}(1.0.9a)

\begin{aligned}T = \left( \frac{\partial {U}}{\partial {S}} \right)_{P} + P \left( \frac{\partial {V}}{\partial {S}} \right)_{P}.\end{aligned} \hspace{\stretch{1}}(1.0.9b)

Differentials in $P, U$

\begin{aligned}0 &= dU - T dS + P dV \\ &= dU - T \left( \left( \frac{\partial {S}}{\partial {P}} \right)_{U} dP + \left( \frac{\partial {S}}{\partial {U}} \right)_{P} dU \right)+ P\left( \left( \frac{\partial {V}}{\partial {P}} \right)_{U} dP + \left( \frac{\partial {V}}{\partial {U}} \right)_{P} dU \right),\end{aligned} \hspace{\stretch{1}}(1.0.10)

or

\begin{aligned}0 = 1 - T \left( \frac{\partial {S}}{\partial {U}} \right)_{P} + P \left( \frac{\partial {V}}{\partial {U}} \right)_{P} \end{aligned} \hspace{\stretch{1}}(1.0.11a)

\begin{aligned}T \left( \frac{\partial {S}}{\partial {P}} \right)_{U} = P \left( \frac{\partial {V}}{\partial {P}} \right)_{U}.\end{aligned} \hspace{\stretch{1}}(1.0.11b)

Differentials in $V, T$

\begin{aligned}0 &= dU - T dS + P dV \\ &= \left( \frac{\partial {U}}{\partial {V}} \right)_{T} dV + \left( \frac{\partial {U}}{\partial {T}} \right)_{V} dT - T \left( \left( \frac{\partial {S}}{\partial {V}} \right)_{T} dV + \left( \frac{\partial {S}}{\partial {T}} \right)_{V} dT \right)+ P dV,\end{aligned} \hspace{\stretch{1}}(1.0.12)

or

\begin{aligned}0 = \left( \frac{\partial {U}}{\partial {V}} \right)_{T} - T \left( \frac{\partial {S}}{\partial {V}} \right)_{T} + P \end{aligned} \hspace{\stretch{1}}(1.0.13a)

\begin{aligned}\left( \frac{\partial {U}}{\partial {T}} \right)_{V} = T \left( \frac{\partial {S}}{\partial {T}} \right)_{V}.\end{aligned} \hspace{\stretch{1}}(1.0.13b)

Differentials in $V, S$

\begin{aligned}0 &= dU - T dS + P dV \\ &= \left( \frac{\partial {U}}{\partial {V}} \right)_{S} dV + \left( \frac{\partial {U}}{\partial {S}} \right)_{V} dS - T dS+ P dV,\end{aligned} \hspace{\stretch{1}}(1.0.14)

or

\begin{aligned}P = -\left( \frac{\partial {U}}{\partial {V}} \right)_{S}\end{aligned} \hspace{\stretch{1}}(1.0.15a)

\begin{aligned}T = \left( \frac{\partial {U}}{\partial {S}} \right)_{V} .\end{aligned} \hspace{\stretch{1}}(1.0.15b)

Differentials in $V, U$

\begin{aligned}0 &= dU - T dS + P dV \\ &= dU- T \left( \left( \frac{\partial {S}}{\partial {V}} \right)_{U} dV + \left( \frac{\partial {S}}{\partial {U}} \right)_{V} dU \right)+ P \left( \left( \frac{\partial {V}}{\partial {V}} \right)_{U} dV + \left( \frac{\partial {V}}{\partial {U}} \right)_{V} dU \right)\end{aligned} \hspace{\stretch{1}}(1.0.16)

or

\begin{aligned}0 = 1 - T \left( \frac{\partial {S}}{\partial {U}} \right)_{V} + P \left( \frac{\partial {V}}{\partial {U}} \right)_{V} \end{aligned} \hspace{\stretch{1}}(1.0.17a)

\begin{aligned}T \left( \frac{\partial {S}}{\partial {V}} \right)_{U} = P \left( \frac{\partial {V}}{\partial {V}} \right)_{U}.\end{aligned} \hspace{\stretch{1}}(1.0.17b)

Differentials in $S, T$

\begin{aligned}0 &= dU - T dS + P dV \\ &= \left( \left( \frac{\partial {U}}{\partial {S}} \right)_{T} dS + \left( \frac{\partial {U}}{\partial {T}} \right)_{S} dT \right)- T dS+ P \left( \left( \frac{\partial {V}}{\partial {S}} \right)_{T} dS + \left( \frac{\partial {V}}{\partial {T}} \right)_{S} dT \right),\end{aligned} \hspace{\stretch{1}}(1.0.18)

or

\begin{aligned}0 = \left( \frac{\partial {U}}{\partial {S}} \right)_{T} - T + P \left( \frac{\partial {V}}{\partial {S}} \right)_{T} \end{aligned} \hspace{\stretch{1}}(1.0.19a)

\begin{aligned}0 = \left( \frac{\partial {U}}{\partial {T}} \right)_{S} + P \left( \frac{\partial {V}}{\partial {T}} \right)_{S}.\end{aligned} \hspace{\stretch{1}}(1.0.19b)

Differentials in $S, U$

\begin{aligned}0 &= dU - T dS + P dV \\ &= dU - T dS+ P \left( \left( \frac{\partial {V}}{\partial {S}} \right)_{U} dS + \left( \frac{\partial {V}}{\partial {U}} \right)_{S} dU \right)\end{aligned} \hspace{\stretch{1}}(1.0.20)

or

\begin{aligned}\frac{1}{{P}} = - \left( \frac{\partial {V}}{\partial {U}} \right)_{S} \end{aligned} \hspace{\stretch{1}}(1.0.21a)

\begin{aligned}T = P \left( \frac{\partial {V}}{\partial {S}} \right)_{U}.\end{aligned} \hspace{\stretch{1}}(1.0.21b)

Differentials in $T, U$

\begin{aligned}0 &= dU - T dS + P dV \\ &= dU - T \left( \left( \frac{\partial {S}}{\partial {T}} \right)_{U} dT + \left( \frac{\partial {S}}{\partial {U}} \right)_{T} dU \right)+ P\left( \left( \frac{\partial {V}}{\partial {T}} \right)_{U} dT + \left( \frac{\partial {V}}{\partial {U}} \right)_{T} dU \right),\end{aligned} \hspace{\stretch{1}}(1.0.22)

or

\begin{aligned}0 = 1 - T \left( \frac{\partial {S}}{\partial {U}} \right)_{T} + P \left( \frac{\partial {V}}{\partial {U}} \right)_{T} \end{aligned} \hspace{\stretch{1}}(1.0.23a)

\begin{aligned}T \left( \frac{\partial {S}}{\partial {T}} \right)_{U} = P \left( \frac{\partial {V}}{\partial {T}} \right)_{U}.\end{aligned} \hspace{\stretch{1}}(1.0.23b)

References

[1] John Baez. Entropic forces, 2012. URL http://johncarlosbaez.wordpress.com/2012/02/01/entropic-forces/. [Online; accessed 07-March-2013].

PHY454H1S Continuum Mechanics. Lecture 15: More on surface tension and Reynold’s number. Taught by Prof. K. Das.

Posted by peeterjoot on March 10, 2012

[Click here for a PDF of this post with nicer formatting and figures if the post had any (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Disclaimer.

Peeter’s lecture notes from class. May not be entirely coherent.

Review. Surface tension clarifications

For a surface like figure (\ref{fig:continuumL15:continuumL15Fig1})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL15Fig1}
\caption{Vapor liquid interface.}
\end{figure}

we have a discontinuous jump in density. We will have to consider three boundary value constraints

\begin{enumerate}
\item Mass balance. This is the continuity equation.
\item Momentum balance. This is the Navier-Stokes equation.
\item Energy balance. This is the heat equation.
\end{enumerate}

We have not yet discussed the heat equation, but this is required for non-isothermal problems.

The boundary condition at the interface is given by the stress balance. Denoting the difference in the traction vector at the interface by

\begin{aligned}[\mathbf{t}]_1^2 = \mathbf{t}_2 - \mathbf{t}_1 = -\frac{\sigma}{2 R} \hat{\mathbf{n}} - \boldsymbol{\nabla}_I \sigma\end{aligned} \hspace{\stretch{1}}(2.1)

Here the gradient is in the tangential direction of the surface as in figure (\ref{fig:continuumL15:continuumL15Fig2})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL15Fig2}
\caption{normal and tangent vectors on a curve.}
\end{figure}

In the normal direction

\begin{aligned}[\mathbf{t}]_1^2 \cdot \hat{\mathbf{n}}&= (\mathbf{t}_2 - \mathbf{t}_1) \cdot \hat{\mathbf{n}} \\ &= -\frac{\sigma}{2 R} \end{aligned}

With the traction vector having the value

\begin{aligned}\mathbf{t} &= \mathbf{e}_i T_{ij} n_j \\ &= \mathbf{e}_i \left( -p \delta_{ij} + \mu \left( \frac{\partial {u_i}}{\partial {x_j}}+\frac{\partial {u_j}}{\partial {x_i}}\right)\right)n_j\end{aligned}

We have in the normal direction

\begin{aligned}\mathbf{t} \cdot \mathbf{n} =n_i \left( -p \delta_{ij} + \mu \left( \frac{\partial {u_i}}{\partial {x_j}}+\frac{\partial {u_j}}{\partial {x_i}}\right)\right) n_j\end{aligned} \hspace{\stretch{1}}(2.2)

With $\mathbf{u} = 0$ on the surface, and $n_i \delta_{ij} n_j = n_j n_j = 1$ we have

\begin{aligned}\mathbf{t} \cdot \mathbf{n} = -p\end{aligned} \hspace{\stretch{1}}(2.3)

Returning to $(\mathbf{t}_2 - \mathbf{t}_1) \cdot \hat{\mathbf{n}}$ we have

\begin{aligned}\boxed{-p_2 + p_1 = -\frac{\sigma}{2 R} }\end{aligned} \hspace{\stretch{1}}(2.4)

This is the Laplace pressure. Note that the sign of the difference is significant, since it effects the direction of the curvature. This is depicted pictorially in figure (\ref{fig:continuumL15:continuumL15Fig3})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL15Fig3}
\caption{pressure and curvature relationships}
\end{figure}

Question In [1] the curvature term is written

\begin{aligned}\frac{1}{{2 R}} \rightarrow \frac{1}{{R_1}} + \frac{1}{{R_2}}.\end{aligned} \hspace{\stretch{1}}(2.5)

Why the difference? Answer: This is to account for non-spherical surfaces, with curvature in two directions. Illustrating by example, imagine a surface like as in figure (\ref{fig:continuumL15:continuumL15Figq})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL15Figq}
\caption{Example of non-spherical curvature.}
\end{figure}

Followup required to truly understand things

While, this review clarifies things, we still don’t really know how the surface tension term $\sigma$ is defined. Nor have we been given any sort of derivation of 2.1, from which the end result follows.

I’m assuming that $\sigma$ is a property of the two fluids at the interface, so if you have, for example, oil and vinegar in a bottle, we have surface tension and curvature that’s probably related to how the two of these interact. If there is still a mixing or settling process occurring, I’d imagine that this could even vary from point to point on the surface (imagine adding soap to a surface where stuff can float until the soap mixes in enough that things start sinking in the radius of influence of the soap).

Surface tension gradients.

Now consider the tangential component of the traction vector

\begin{aligned}\mathbf{t}_2 \cdot \hat{\boldsymbol{\tau}} - \mathbf{t}_1 \cdot \hat{\boldsymbol{\tau}} = - \not{{ \frac{\sigma}{2 R} \hat{\mathbf{n}} \cdot \hat{\boldsymbol{\tau}}}} - \hat{\boldsymbol{\tau}} \cdot \boldsymbol{\nabla}_I \sigma\end{aligned} \hspace{\stretch{1}}(2.6)

So we see that for a static fluid, we must have

\begin{aligned}\boldsymbol{\nabla}_I \sigma = 0\end{aligned} \hspace{\stretch{1}}(2.7)

For a static interface there cannot be any surface tension gradient. This becomes very important when considering stability issues. We can have surface tension induced flow called capillary, or mandarin (?) flow.

Reynold’s number.

In Navier-Stokes after making non-dimensionalization changes of the form

\begin{aligned}x \rightarrow L x'\end{aligned} \hspace{\stretch{1}}(3.8)

the control parameter is like Reynold’s number.

In NS

\begin{aligned}\rho \frac{\partial {\mathbf{u}}}{\partial {t}} + \rho ( \mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u} = - \boldsymbol{\nabla} p + \mu \boldsymbol{\nabla}^2 \mathbf{u}\end{aligned} \hspace{\stretch{1}}(3.9)

We call the term

\begin{aligned}\rho ( \mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u}\end{aligned} \hspace{\stretch{1}}(3.10)

the inertial term. It is non-zero only when something is being “carried along with the velocity”. Consider a volume fixed in space and one that is moving along with the fluid as in figure (\ref{fig:continuumL15:continuumL15Fig4})
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL15Fig4}
\caption{Moving and fixed frame control volumes in a fluid.}
\end{figure}

All of our viscosity dependence shows up in the Laplacian term, so we can roughly characterize the Reynold’s number as the ratio

\begin{aligned}\text{Reynold's number}&\rightarrow\frac{{\left\lvert{\text{effect of inertia}}\right\rvert}}{{\left\lvert{\text{effect of viscosity}}\right\rvert}} \\ &=\frac{ {\left\lvert{\rho ( \mathbf{u} \cdot \boldsymbol{\nabla} ) \mathbf{u} }\right\rvert} }{{\left\lvert{ \mu \boldsymbol{\nabla}^2 \mathbf{u}}\right\rvert}} \\ &\sim\frac{ \rho U^2/L }{\mu U/L^2} \\ &\sim\frac{ \rho U L }{\mu }\end{aligned}

In figures (\ref{fig:continuumL15:continuumL15Fig5}), (\ref{fig:continuumL15:continuumL15Fig6}) we have two illustrations of viscous and non-viscous regions the first with a moving probe pushing its way through a surface, and the second with a wing set at an angle of attack that generates some turbulence. Both are illustrations of the viscous and inviscous regions for the two flows. Both of these are characterized by the Reynold’s number in some way not really specified in class. One of the points of mentioning this is that when we are in an essentially inviscous region, we can neglect the viscosity ($\mu \boldsymbol{\nabla}^2 \mathbf{u}$) term of the flow.
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL15Fig5}
\caption{continuumL15Fig5}
\end{figure}
\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{continuumL15Fig6}
\caption{continuumL15Fig6}
\end{figure}

References

[1] L.D. Landau and E.M. Lifshitz. A Course in Theoretical Physics-Fluid Mechanics. Pergamon Press Ltd., 1987.