Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Posts Tagged ‘Gamma function’

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.

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , | 1 Comment »

Non integral binomial coefficient

Posted by peeterjoot on April 10, 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.)]

Updated.  Original generalization of binomial coefficient wasn’t correct for negative exponents.

In [2] appendix section F was the use of binomial coefficients in a non-integral binomial expansion. This surprised me, since I’d never seen that before. However, on reflection, this is a very sensible notation, provided the binomial coefficients are defined in terms of the gamma function. Let’s explore this little detail explicitly.

Taylor series

We start with a Taylor expansion of

\begin{aligned}f(x) = (a + x)^b.\end{aligned} \hspace{\stretch{1}}(1.1)

Our derivatives are

\begin{aligned}\begin{aligned}f'(x) &= b (a + x)^{b-1} \\ f''(x) &= b(b-1) (a + x)^{b-2} \\ f^3(x) &= b(b-1)(b-(3-1)) (a + x)^{b-3} \\ \dot{v}s & \\ f^k(x) &= b(b-1)\cdots (b-(k-1)) (a + x)^{b-k}.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.2)

Our Taylor series is then

\begin{aligned}(a + x)^b = \sum_{k = 0}^\infty \frac{1}{{k!}} b(b-1)\cdots (b-(k-1)) a^{b-k} x^k.\end{aligned} \hspace{\stretch{1}}(1.0.2)

Note that if b is a positive integer, then all the elements of this series become zero at b = k - 1, or

\begin{aligned}(a + x)^b = \sum_{k = 0}^k \frac{1}{{k!}} b(b-1)\cdots (b-(k-1)) a^{b-k} x^k.\end{aligned} \hspace{\stretch{1}}(1.0.2)

Gamma function

Let’s now relate this to the gamma function. From [1] section 6.1.1 we have

\begin{aligned}\Gamma(z) = \int_0^\infty t^{z-1} e^{-t} dt.\end{aligned} \hspace{\stretch{1}}(1.0.2)

Iteratively integrating by parts, we find the usual relation between gamma functions of integral separation

\begin{aligned}\Gamma(z + 1) &= \int_0^\infty t^{z} e^{-t} dt \\ &= \int_0^\infty t^{z} d \left( \frac{ e^{-t}}{-1} \right) \\ &= {\left.{{t^z \frac{e^{-t}}{-1}}}\right\vert}_{{0}}^{{\infty}}- \int_0^\infty z t^{z-1} \frac{e^{-t}}{-1} dt \\ &= z \int_0^\infty t^{z-1} e^{-t} dt \\ &= z (z - 1)\int_0^\infty t^{z-2} e^{-t} dt \\ &= z (z - 1)(z - (3-1))\int_0^\infty t^{z-3} e^{-t} dt \\ &= z (z - 1) \cdots (z - (k-1))\int_0^\infty t^{z-k} e^{-t} dt \\ &= z (z - 1) \cdots (z - (k-1))\int_0^\infty t^{(z + 1 - k) - 1} e^{-t} dt,\end{aligned} \hspace{\stretch{1}}(1.0.2)


\begin{aligned}\Gamma(z + 1) = z (z - 1) \cdots (z - (k-1)) \Gamma( z - (k-1) ).\end{aligned} \hspace{\stretch{1}}(1.0.2)

Flipping this gives us a nice closed form expression for the products of a number of positive unit separated values

\begin{aligned}z (z - 1) \cdots (z - k)=\frac{\Gamma(z + 1)}{\Gamma( z - (k-1) )}.\end{aligned} \hspace{\stretch{1}}(1.0.2)

Binomial coefficient for positive exponents

Considering first positive exponents b, we can now use this in our Taylor expansion eq. 1.0.2

\begin{aligned}(a + x)^b &= \sum_{k = 0}^\infty \frac{1}{{k!}} \frac{\Gamma(b+1)}{\Gamma(b - k + 1)}a^{b-k} x^k \\ &= \sum_{k = 0}^\infty \frac{\Gamma(b + 1)}{\Gamma(k + 1)\Gamma(b - k + 1)}a^{b-k} x^k.\end{aligned} \hspace{\stretch{1}}(1.0.2)

Observe that when b is a positive integer we have

\begin{aligned}\frac{\Gamma(b + 1)}{\Gamma(k + 1)\Gamma(b - k + 1)} &= \frac{b!}{k!(b-k)!} \\ &= \binom{b}{k}.\end{aligned} \hspace{\stretch{1}}(1.0.10)

So for positive values of b, even non-integer values, we see that is then very reasonable to define the binomial coefficient \index{binomial coefficient} explicitly in terms of the gamma function

\begin{aligned}\binom{b}{k} \equiv\frac{\Gamma(b + 1)}{\Gamma(k + 1)\Gamma(b - k + 1)}.\end{aligned} \hspace{\stretch{1}}(1.0.11)

If we do that, then the binomial expansion for non-integral values of b is simply

\begin{aligned}(a + x)^b = \sum_{k = 0}^\infty \binom{b}{k} a^{b-k} x^k.\end{aligned} \hspace{\stretch{1}}(1.0.11)

Binomial coefficient for negative integer exponents

Using the relation eq. 1.0.11 blindly leads to some trouble, since \Gamma(-\left\lvert {m} \right\rvert) goes to infinity for integer values of m > 0. We have to modify the definition of the binomial coefficient. Let’s rewrite eq. 1.0.2 for negative integer values of b = -m as

\begin{aligned}(a + x)^{-m} &= \sum_{k = 0}^\infty \frac{1}{{k!}} (-m)(-m-1)\cdots (-m-(k-1)) a^{-m-k} x^k \\ &= \sum_{k = 0}^\infty \frac{1}{{k!}} (-1)^km(m+1)\cdots (m+(k-1)) a^{-m-k} x^k.\end{aligned} \hspace{\stretch{1}}(1.0.11)

Let’s also put the ratio of gamma functions relation of eq. 1.0.2, in a slightly more general form. For u, v > 0, where u - v is an integer, we can write

\begin{aligned}u (u - 1) \cdots (v) = \frac{\Gamma(u + 1)}{\Gamma( v )}.\end{aligned} \hspace{\stretch{1}}(1.0.11)

Our Taylor series takes the form

\begin{aligned}(a + x)^{-m} &= \sum_{k = 0}^\infty \frac{1}{{k!}} (-1)^k(m+k-1) (m+k-2) \cdots (m)a^{-m-k} x^k \\ &= \sum_{k = 0}^\infty (-1)^k \frac{ \Gamma(m + k) }{\Gamma(m)\Gamma(k+1)}a^{-m-k} x^k.\end{aligned} \hspace{\stretch{1}}(1.0.16)

We can now define, for negative integers -m

\begin{aligned}\binom{-m}{k}\equiv(-1)^k \frac{ \Gamma(m + k) }{ \Gamma(m)\Gamma(k+1) }.\end{aligned} \hspace{\stretch{1}}(1.0.16)

With such a definition, our Taylor series takes the tidy form

\begin{aligned}(a + x)^{-m} = \sum_{k = 0}^\infty \binom{-m}{k} a^{-m-k} x^k.\end{aligned} \hspace{\stretch{1}}(1.0.16)

For negative integer values of b = -m, this is now consistent with eq. 1.0.11.

Observe that we can put eq. 1.0.16 into the standard binomial form with a bit of manipulation

\begin{aligned}\binom{-m}{k} &= (-1)^k \frac{ \Gamma(m + k) }{ \Gamma(m)\Gamma(k+1) } \\ &= (-1)^k \frac{ (m + k -1)! }{ (m-1)! k! } \\ &= (-1)^k \frac{ m (m + k)! }{ (m + k) m! k! },\end{aligned} \hspace{\stretch{1}}(1.0.16)


\begin{aligned}\binom{-m}{k}=(-1)^k \frac{m}{m + k} \binom{m+k}{m}.\end{aligned} \hspace{\stretch{1}}(1.0.19)

Negative non-integral binomial coefficients

TODO. There will be some ugliness due to the changes of sign in the products b(b-1)\cdots (b -k + 1) since b and b -k + 1 may have different sign. A product of two ratios of gamma functions will be required to express this product, which will further complicate the definition of binomial coefficient.


[1] M. Abramowitz and I.A. Stegun. \emph{Handbook of mathematical functions with formulas, graphs, and mathematical tables}, volume 55. Dover publications, 1964.

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

Posted in Math and Physics Learning. | Tagged: , , , , | Leave a Comment »

PHY452H1S Basic Statistical Mechanics. Problem Set 6: Max entropy, fugacity, and Fermi gas

Posted by peeterjoot on April 3, 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.)]


This is an ungraded set of answers to the problems posed.

Question: Maximum entropy principle

Consider the “Gibbs entropy”

\begin{aligned}S = - k_{\mathrm{B}} \sum_i p_i \ln p_i\end{aligned} \hspace{\stretch{1}}(1.1)

where p_i is the equilibrium probability of occurrence of a microstate i in the ensemble.

For a microcanonical ensemble with \Omega configurations (each having the same energy), assigning an equal probability p_i= 1/\Omega to each microstate leads to S = k_{\mathrm{B}} \ln \Omega. Show that this result follows from maximizing the Gibbs entropy with respect to the parameters p_i subject to the constraint of

\begin{aligned}\sum_i p_i = 1\end{aligned} \hspace{\stretch{1}}(1.2)

(for p_i to be meaningful as probabilities). In order to do the minimization with this constraint, use the method of Lagrange multipliers – first, do an unconstrained minimization of the function

\begin{aligned}S - \alpha \sum_i p_i,\end{aligned} \hspace{\stretch{1}}(1.3)

then fix \alpha by demanding that the constraint be satisfied.

For a canonical ensemble (no constraint on total energy, but all microstates having the same number of particles N), maximize the Gibbs entropy with respect to the parameters p_i subject to the constraint of

\begin{aligned}\sum_i p_i = 1,\end{aligned} \hspace{\stretch{1}}(1.4)

(for p_i to be meaningful as probabilities) and with a given fixed average energy

\begin{aligned}\left\langle{{E}}\right\rangle = \sum_i E_i p_i,\end{aligned} \hspace{\stretch{1}}(1.5)

where E_i is the energy of microstate i. Use the method of Lagrange multipliers, doing an unconstrained minimization of the function

\begin{aligned}S - \alpha \sum_i p_i - \beta \sum_i E_i p_i,\end{aligned} \hspace{\stretch{1}}(1.6)

then fix \alpha, \beta by demanding that the constraint be satisfied. What is the resulting p_i?

For a grand canonical ensemble (no constraint on total energy, or the number of particles), maximize the Gibbs entropy with respect to the parameters p_i subject to the constraint of

\begin{aligned}\sum_i p_i = 1,\end{aligned} \hspace{\stretch{1}}(1.7)

(for p_i to be meaningful as probabilities) and with a given fixed average energy

\begin{aligned}\left\langle{{E}}\right\rangle = \sum_i E_i p_i,\end{aligned} \hspace{\stretch{1}}(1.8)

and a given fixed average particle number

\begin{aligned}\left\langle{{N}}\right\rangle = \sum_i N_i p_i.\end{aligned} \hspace{\stretch{1}}(1.9)

Here E_i, N_i represent the energy and number of particles in microstate i. Use the method of Lagrange multipliers, doing an unconstrained minimization of the function

\begin{aligned}S - \alpha \sum_i p_i - \beta \sum_i E_i p_i - \gamma \sum_i N_i p_i,\end{aligned} \hspace{\stretch{1}}(1.10)

then fix \alpha, \beta, \gamma by demanding that the constrains be satisfied. What is the resulting p_i?



\begin{aligned}f = S - \alpha \sum_{j = 1}^\Omega p_j,= -\sum_{j = 1}^\Omega p_j \left( k_{\mathrm{B}} \ln p_j + \alpha  \right),\end{aligned} \hspace{\stretch{1}}(1.11)

our unconstrained minimization requires

\begin{aligned}0 = \frac{\partial {f}}{\partial {p_i}}= -\left( k_{\mathrm{B}} \left( \ln p_i + 1  \right) + \alpha  \right).\end{aligned} \hspace{\stretch{1}}(1.12)

Solving for p_i we have

\begin{aligned}p_i = e^{-\alpha/k_{\mathrm{B}} - 1}.\end{aligned} \hspace{\stretch{1}}(1.13)

The probabilities for each state are constant. To fix that constant we employ our constraint

\begin{aligned}1 = \sum_{j = 1}^\Omega p_j= \sum_{j = 1}^\Omega e^{-\alpha/k_{\mathrm{B}} - 1}= \Omega e^{-\alpha/k_{\mathrm{B}} - 1},\end{aligned} \hspace{\stretch{1}}(1.14)


\begin{aligned}\alpha/k_{\mathrm{B}} + 1 = \ln \Omega.\end{aligned} \hspace{\stretch{1}}(1.15)

Inserting eq. 1.15 fixes the probability, giving us the first of the expected results

\begin{aligned}\boxed{p_i = e^{-\ln \Omega} = \frac{1}{{\Omega}}.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

Using this we our Gibbs entropy can be summed easily

\begin{aligned}S &= -k_{\mathrm{B}} \sum_{j = 1}^\Omega p_j \ln p_j \\ &= -k_{\mathrm{B}} \sum_{j = 1}^\Omega \frac{1}{{\Omega}} \ln \frac{1}{{\Omega}} \\ &= -k_{\mathrm{B}} \frac{\Omega}{\Omega} \left( -\ln \Omega  \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)


\begin{aligned}\boxed{S = k_{\mathrm{B}} \ln \Omega.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

For the “action” like quantity that we want to minimize, let’s write

\begin{aligned}f = S - \alpha \sum_j p_j - \beta \sum_j E_j p_j,\end{aligned} \hspace{\stretch{1}}(1.0.16)

for which we seek \alpha, \beta such that

\begin{aligned}0 &= \frac{\partial {f}}{\partial {p_i}} \\ &= -\frac{\partial {}}{\partial {p_i}}\sum_j p_j\left( k_{\mathrm{B}} \ln p_j + \alpha + \beta E_j  \right) \\ &= -k_{\mathrm{B}} (\ln p_i + 1) - \alpha - \beta E_i,\end{aligned} \hspace{\stretch{1}}(1.0.16)


\begin{aligned}p_i = \exp\left( - \left( \alpha - \beta E_i \right) /k_{\mathrm{B}} - 1  \right).\end{aligned} \hspace{\stretch{1}}(1.0.16)

Our probability constraint is

\begin{aligned}1 &= \sum_j \exp\left( - \left( \alpha - \beta E_j \right) /k_{\mathrm{B}} - 1  \right) \\ &= \exp\left( - \alpha/k_{\mathrm{B}} - 1  \right)\sum_j \exp\left( - \beta E_j/k_{\mathrm{B}}  \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)


\begin{aligned}\exp\left( \alpha/k_{\mathrm{B}} + 1  \right)=\sum_j \exp\left( - \beta E_j/k_{\mathrm{B}}  \right).\end{aligned} \hspace{\stretch{1}}(1.0.16)

Taking logs we have

\begin{aligned}\alpha/k_{\mathrm{B}} + 1 = \ln \sum_j \exp\left( - \beta E_j/k_{\mathrm{B}}  \right).\end{aligned} \hspace{\stretch{1}}(1.0.16)

We could continue to solve for \alpha explicitly but don’t care any more than this. Plugging back into the probability eq. 1.0.16 obtained from the unconstrained minimization we have

\begin{aligned}p_i = \exp\left( -\ln \sum_j \exp\left( - \beta E_j/k_{\mathrm{B}}  \right)  \right)\exp\left( - \beta E_i/k_{\mathrm{B}}  \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)


\begin{aligned}\boxed{p_i = \frac{   \exp\left( - \beta E_i/k_{\mathrm{B}}  \right)}{   \sum_j    \exp\left( - \beta E_j/k_{\mathrm{B}}  \right)}.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

To determine \beta we must look implicitly to the energy constraint, which is

\begin{aligned}\left\langle{{E}}\right\rangle &= \sum_i E_i p_i \\ &= \sum_iE_i\left( \frac{ \exp\left( - \beta E_i/k_{\mathrm{B}}  \right) } { \sum_j \exp\left( - \beta E_j/k_{\mathrm{B}}  \right) }  \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)


\begin{aligned}\boxed{\left\langle{{E}}\right\rangle = \frac{   \sum_i E_i \exp\left( - \beta E_i/k_{\mathrm{B}}  \right)}{   \sum_j \exp\left( - \beta E_j/k_{\mathrm{B}}  \right)}.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

The constraint \beta (=1/T) is given implicitly by this energy constraint.

Again write

\begin{aligned}f = S - \alpha \sum_j p_j - \beta \sum_j E_j p_j - \gamma \sum_j N_j p_j.\end{aligned} \hspace{\stretch{1}}(1.0.16)

The unconstrained minimization requires

\begin{aligned}0 &= \frac{\partial {f}}{\partial {p_i}} \\ &= -\frac{\partial {}}{\partial {p_i}}\left( k_{\mathrm{B}} (\ln p_i + 1) + \alpha + \beta E_i + \gamma N_i  \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)


\begin{aligned}p_i = \exp\left( -\alpha/k_{\mathrm{B}} - 1  \right) \exp\left( -(\beta E_i + \gamma N_i)/k_{\mathrm{B}}  \right).\end{aligned} \hspace{\stretch{1}}(1.0.16)

The unit probability constraint requires

\begin{aligned}1 &= \sum_j p_j \\ &= \exp\left( -\alpha/k_{\mathrm{B}} - 1  \right) \sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}}  \right),\end{aligned} \hspace{\stretch{1}}(1.0.16)


\begin{aligned}\exp\left( -\alpha/k_{\mathrm{B}} - 1  \right) =\frac{1}{{\sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}}  \right)}}.\end{aligned} \hspace{\stretch{1}}(1.0.16)

Our probability is then

\begin{aligned}\boxed{p_i = \frac{\exp\left( -(\beta E_i + \gamma N_i)/k_{\mathrm{B}}  \right)}{\sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}}  \right)}.}\end{aligned} \hspace{\stretch{1}}(1.0.16)

The average energy \left\langle{{E}}\right\rangle = \sum_j p_j E_j and average number of particles \left\langle{{N}}\right\rangle = \sum_j p_j N_j are given by

\begin{aligned}\left\langle{{E}}\right\rangle = \frac{E_i\exp\left( -(\beta E_i + \gamma N_i)/k_{\mathrm{B}}  \right)}{\sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}}  \right)}\end{aligned} \hspace{\stretch{1}}(

\begin{aligned}\left\langle{{N}}\right\rangle = \frac{N_i\exp\left( -(\beta E_i + \gamma N_i)/k_{\mathrm{B}}  \right)}{\sum_j\exp\left( -(\beta E_j + \gamma N_j)/k_{\mathrm{B}}  \right)}.\end{aligned} \hspace{\stretch{1}}(

The values \beta and \gamma are fixed implicitly by requiring simultaneous solutions of these equations.

Question: Fugacity expansion ([3] Pathria, Appendix D, E)

The theory of the ideal Fermi or Bose gases often involves integrals of the form

\begin{aligned}f_\nu^\pm(z) = \frac{1}{{\Gamma(\nu)}} \int_0^\infty dx \frac{x^{\nu - 1}}{z^{-1} e^x \pm 1}\end{aligned} \hspace{\stretch{1}}(1.36)


\begin{aligned}\Gamma(\nu) = \int_0^\infty dy y^{\nu-1} e^{-y}\end{aligned} \hspace{\stretch{1}}(1.37)

denotes the gamma function.

Obtain the behavior of f_\nu^\pm(z) for z \rightarrow 0 keeping the two leading terms in the expansion.

For Fermions, obtain the behavior of f_\nu^\pm(z) for z \rightarrow \infty again keeping the two leading terms.

For Bosons, we must have z \le 1 (why?), obtain the leading term of f_\nu^-(z) for z \rightarrow 1.


For z \rightarrow 0 we can rewrite the integrand in a form that allows for series expansion

\begin{aligned}\frac{x^{\nu - 1}}{z^{-1} e^x \pm 1} &= \frac{z e^{-x} x^{\nu - 1}}{1 \pm z e^{-x}} \\ &= z e^{-x} x^{\nu - 1}\left( 1 \mp z e^{-x} + (z e^{-x})^2 \mp (z e^{-x})^3 + \cdots  \right)\end{aligned} \hspace{\stretch{1}}(1.38)

For the kth power of z e^{-x} in this series our integral is

\begin{aligned}\int_0^\infty dx z e^{-x} x^{\nu - 1} (z e^{-x})^k &= z^{k+1}\int_0^\infty dx x^{\nu - 1} e^{-(k + 1) x} \\ &= \frac{z^{k+1}}{(k+1)^\nu}\int_0^\infty du u^{\nu - 1} e^{- u} \\ &= \frac{z^{k+1}}{(k+1)^\nu} \Gamma(\nu)\end{aligned} \hspace{\stretch{1}}(1.39)

Putting everything back together we have for small z

\begin{aligned}\boxed{f_\nu^\pm(z) =z\mp\frac{z^{2}}{2^\nu}+\frac{z^{3}}{3^\nu}\mp\frac{z^{4}}{4^\nu}+ \cdots}\end{aligned} \hspace{\stretch{1}}(1.40)

We’ll expand \Gamma(\nu) f_\nu^+(e^y) about z = e^y, writing

\begin{aligned}F_\nu(e^y) &= \Gamma(\nu) f_\nu^+(e^y) \\ &= \int_0^\infty dx \frac{x^{\nu - 1}}{e^{x - y} + 1} \\ &= \int_0^y dx \frac{x^{\nu - 1}}{e^{x - y} + 1}+\int_y^\infty dx \frac{x^{\nu - 1}}{e^{x - y} + 1}.\end{aligned} \hspace{\stretch{1}}(1.41)

The integral has been split into two since the behavior of the exponential in the denominator is quite different in the x  y ranges. Observe that in the first integral we have

\begin{aligned}\frac{1}{{2}} \le \frac{1}{e^{x - y} + 1} \le \frac{1}{{1 + e^{-y}}}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

Since this term is of order 1, let’s consider the difference of this from 1, writing

\begin{aligned}\frac{1}{e^{x - y} + 1} = 1 + u,\end{aligned} \hspace{\stretch{1}}(1.0.42)


\begin{aligned}u = \frac{1}{e^{x - y} + 1} - 1 &= \frac{1 -(e^{x - y} + 1)}{e^{x - y} + 1} \\ &= \frac{-e^{x - y} }{e^{x - y} + 1} \\ &= -\frac{1}{1 + e^{y - x}}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

This gives us

\begin{aligned}F_\nu(e^y) &= \int_0^y dx x^{\nu - 1} \left( 1 - \frac{ 1 } { 1 + e^{y - x} }  \right)+\int_y^\infty dx \frac{x^{\nu - 1}}{e^{x - y} + 1} \\ &= \frac{y^\nu}{\nu}-\int_0^y dx \frac{ x^{\nu - 1}  } { 1 + e^{y - x} }+\int_y^\infty dx \frac{x^{\nu - 1}}{e^{x - y} + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

Now let’s make a change of variables a = y - x in the first integral and b = x - y in the second. This gives

\begin{aligned}F_\nu(e^y) = \frac{y^\nu}{\nu}-\int_0^\infty da \frac{ (y - a)^{\nu - 1}  } { 1 + e^{a} }+\int_0^\infty db \frac{(y + b)^{\nu - 1}}{e^{b} + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

As a gets large in the first integral the integrand is approximately e^{-a} (y-a)^{\nu -1}. The exponential dominates this integrand. Since we are considering large y, we can approximate the upper bound of the integral by extending it to \infty. Also expanding in series we have

\begin{aligned}F_\nu(e^y)  &\approx \frac{y^\nu}{\nu}+\int_0^\infty da \frac{ (y + a)^{\nu - 1} -(y - a)^{\nu - 1}  } { 1 + e^{a} } \\ &= \frac{y^\nu}{\nu}+\int_0^\infty da \frac{1}{{e^a + 1}}\left( \left( \frac{1}{{0!}} y^{\nu-1} a^0 + \frac{1}{{1!}} (\nu-1) y^{\nu-2} a^1 + \frac{1}{{2!}} (\nu-1) (\nu-2) y^{\nu-3} a^2 + \cdots  \right) - \left( \frac{1}{{0!}} y^{\nu-1} (-a)^0 + \frac{1}{{1!}} (\nu-1) y^{\nu-2} (-a)^1 + \frac{1}{{2!}} (\nu-1) (\nu-2) y^{\nu-3} (-a)^2 + \cdots  \right)  \right) \\ &= \frac{y^\nu}{\nu}+ 2\int_0^\infty da \frac{1}{{e^a + 1}}   \left( \frac{1}{{1!}} (\nu-1) y^{\nu-2} a^1 + \frac{1}{{3!}} (\nu-1) (\nu-2) (\nu - 3)y^{\nu-4} a^3 + \cdots  \right) \\ &= \frac{y^\nu}{\nu}+ 2\sum_{j = 1, 3, 5, \cdots} \frac{y^{\nu - 1 - j}}{j!} \left( \prod_{k = 1}^j (\nu-k)  \right)\int_0^\infty da \frac{a^j}{e^a + 1} \\ &= \frac{y^\nu}{\nu}+ 2\sum_{j = 1, 3, 5, \cdots} \frac{y^{\nu - 1 - j}}{j!} \frac{ \Gamma(\nu) } {\Gamma(\nu - j)}\int_0^\infty da \frac{a^j}{e^a + 1}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

For the remaining integral, Mathematica gives

\begin{aligned}\int_0^\infty da \frac{a^j}{e^a + 1}=\left( 1 - 2^{-j} \right) j! \zeta (j+1),\end{aligned} \hspace{\stretch{1}}(1.0.42)

where for s > 1

\begin{aligned}\zeta(s) = \sum_{k=1}^{\infty} k^{-s}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

This gives

\begin{aligned}F_\nu(e^y)  \approx \frac{y^\nu}{\nu}+ 2\sum_{j = 1, 3, 5, \cdots} y^{\nu - 1 - j}\frac{ \Gamma(\nu) } {\Gamma(\nu - j)}\left( 1 - 2^{-j} \right) \zeta(j+1),\end{aligned} \hspace{\stretch{1}}(1.0.42)


\begin{aligned}f_\nu^+(e^y)  &\approx y^\nu\left( \frac{1}{\nu \Gamma(\nu)} + 2 \sum_{j = 1, 3, 5, \cdots} \frac{ 1 } {\Gamma(\nu - j) y^{j + 1} } \left( 1 - 2^{-j} \right) \zeta(j+1)  \right) \\ &= \frac{y^\nu}{\Gamma(\nu + 1)}\left( 1 + 2 \sum_{j = 1, 3, 5, \cdots} \frac{ \Gamma(\nu + 1) } {\Gamma(\nu - j) } \left( 1 - 2^{-j} \right) \frac{\zeta(j+1)}{ y^{j + 1} }  \right),\end{aligned} \hspace{\stretch{1}}(1.0.42)


\begin{aligned}\boxed{f_\nu^+(e^y)  \approx \frac{y^\nu}{\Gamma(\nu + 1)}\left( 1 + 2 \nu \sum_{j = 1, 3, 5, \cdots} (\nu-1) \cdots(\nu - j) \left( 1 - 2^{-j} \right) \frac{\zeta(j+1)}{ y^{j + 1} }  \right).}\end{aligned} \hspace{\stretch{1}}(1.0.42)

Evaluating the numerical portions explicitly, with

\begin{aligned}c(j) = 2 \left(1-2^{-j}\right) \zeta (j+1),\end{aligned} \hspace{\stretch{1}}(1.0.42)

\begin{aligned}\begin{aligned}c(1) &= \frac{\pi^2}{6} \\ c(3) &= \frac{7 \pi^4}{360} \\ c(5) &= \frac{31 \pi^6}{15120} \\ c(7) &= \frac{127 \pi^8}{604800},\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.42)

so to two terms (j = 1, 3), we have

\begin{aligned}\boxed{f_\nu^+(e^y)  \approx \frac{y^\nu}{\Gamma(\nu + 1)}\left( 1 + \nu(\nu-1) \frac{\pi^2}{6 y^{2}} + \nu(\nu-1)(\nu-2)(\nu -3) \frac{7 \pi^4}{360 y^4}  \right).}\end{aligned} \hspace{\stretch{1}}(1.0.42)

In order for the Boson occupation numbers to be non-singular we require \mu less than all \epsilon. If that lowest energy level is set to zero, this is equivalent to z < 1. Given this restriction, a z = e^{-\alpha} substitution is convenient for investigation of the z \rightarrow 1 case. Following the text, we'll write

\begin{aligned}G_\nu(e^{-\alpha})=\Gamma(\nu)f_\nu^-(e^{-\alpha}) = \int_0^\infty dx \frac{x^{\nu - 1}}{e^{x + \alpha} - 1}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

For \nu = 1, this is integrable

\begin{aligned}\frac{d}{dx} \ln\left( 1 - e^{-x - \alpha}  \right) &= \frac{e^{-x - \alpha}}{ 1 - e^{-x - \alpha} } \\ &= \frac{1}{{ e^{x + \alpha} - 1}},\end{aligned} \hspace{\stretch{1}}(1.0.42)

so that

\begin{aligned}G_1(e^{-\alpha}) &= \int_0^\infty dx \frac{1}{e^{x + \alpha} - 1} \\ &= {\ln \left( 1 - e^{-x - \alpha}  \right)}_{0}^{\infty} \\ &= \ln 1 - \ln \left( 1 - e^{- \alpha}  \right) \\ &= -\ln \left( 1 - e^{- \alpha}  \right).\end{aligned} \hspace{\stretch{1}}(1.0.42)

Taylor expanding 1 - e^{-\alpha} we have

\begin{aligned}1 - e^{-\alpha} = 1 - \left( 1 - \alpha + \alpha^2/2 - \cdots \right).\end{aligned} \hspace{\stretch{1}}(1.0.42)

Noting that \Gamma(1) = 1, we have for the limit

\begin{aligned}\lim_{\alpha \rightarrow 0} G_1(e^{-\alpha}) \rightarrow - \ln \alpha,\end{aligned} \hspace{\stretch{1}}(1.0.42)


\begin{aligned}\lim_{z\rightarrow 1} f_\nu^-(z)= -\ln (-\ln z).\end{aligned} \hspace{\stretch{1}}(1.0.42)

For values of \nu \ne 1, the denominator is

\begin{aligned}e^{\alpha + x} - 1 = (\alpha + x) + (\alpha + x)^2/2 + \cdots\end{aligned} \hspace{\stretch{1}}(1.0.42)

To first order this gives us

\begin{aligned}f_\nu^-( e^{-\alpha} ) \approx \frac{1}{{\Gamma(\nu)}} \int_0^\infty dx \frac{1}{x + \alpha}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

Of this integral Mathematica says it can be evaluated for 0 < \nu < 1, and has the value

\begin{aligned}\frac{1}{{\Gamma(\nu)}} \int_0^\infty dx \frac{1}{x + \alpha}=\frac{\pi}{\sin(\pi\nu)} \frac{1}{\alpha^{1 - \nu} \Gamma (\nu )}.\end{aligned} \hspace{\stretch{1}}(1.0.42)

From [1] 6.1.17 we find

\begin{aligned}\Gamma(z) \Gamma(1-z) = \frac{\pi}{\sin(\pi z)},\end{aligned} \hspace{\stretch{1}}(1.0.42)

with which we can write

\begin{aligned}\boxed{f_\nu^-( e^{-\alpha} ) \approx \frac{ \Gamma(1 - \nu)}{ \alpha^{1 - \nu} }.}\end{aligned} \hspace{\stretch{1}}(1.0.42)

Question: Nuclear matter ([2], prob 9.2)

Consider a heavy nucleus of mass number A. i.e., having A total nucleons including neutrons and protons. Assume that the number of neutrons and protons is equal, and recall that each of them has spin-1/2 (so possessing two spin states). Treating these nucleons as a free ideal Fermi gas of uniform density contained in a radius R = r_0 A^{1/3}, where r_0 = 1.4 \times 10^{-13} \text{cm}, calculate the Fermi energy and the average energy per nucleon in MeV.


Our nucleon particle density is

\begin{aligned}\rho &= \frac{N}{V} \\ &= \frac{A}{\frac{4 \pi}{3} R^3} \\ &= \frac{3 A}{4 \pi r_0^3 A} \\ &= \frac{3}{4 \pi r_0^3} \\ &= \frac{3}{4 \pi (1.4 \times 10^{-13} \text{cm})^3} \\ &= 8.7 \times 10^{37} (\text{cm})^{-3} \\ &= 8.7 \times 10^{43} (\text{m})^{-3}\end{aligned} \hspace{\stretch{1}}(1.67)

With m for the mass of either the proton or the neutron, and \rho_m = \rho_p = \rho/2, the Fermi energy for these particles is

\begin{aligned}\epsilon_{\mathrm{F}} = \frac{\hbar^2}{2m} \left( \frac{6 \pi (\rho/2)}{2 S + 1} \right)^{2/3},\end{aligned} \hspace{\stretch{1}}(1.68)

With S = 1/2, and 2 S + 1 = 2(1/2) + 1 = 2 for either the proton or the neutron, this is

\begin{aligned}\epsilon_{\mathrm{F}} = \frac{\hbar^2}{2 m} \left( \frac{3 \pi^2 \rho}{2} \right)^{2/3}.\end{aligned} \hspace{\stretch{1}}(1.69)

\begin{aligned}\begin{aligned}\hbar &= 1.05 \times 10^{-34} \,\text{m^2 kg s^{-1}} \\ m &= 1.67 \times 10^{-27} \,\text{kg}\end{aligned}.\end{aligned} \hspace{\stretch{1}}(1.70)

This gives us

\begin{aligned}\epsilon_{\mathrm{F}} &= \frac{(1.05 \times 10^{-34})^2}{2 \times 1.67 \times 10^{-27}} \left( \frac{3 \pi^2 }{2} \frac{8.7 \times 10^{43} }{2}  \right)^{2/3}\text{m}^4 \frac{\text{kg}^2}{s^2} \frac{1}{{\text{kg}}} \frac{1}{{\text{m}^2}} \\ &= 3.9 \times 10^{-12} \,\text{J} \times \left( 6.241509 \times 10^{12} \frac{\text{MeV}}{J} \right) \approx 24 \text{MeV}\end{aligned} \hspace{\stretch{1}}(1.71)

In lecture 16

we found that the total average energy for a Fermion gas of N particles was

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

so the average energy per nucleon is approximately

\begin{aligned}\frac{3}{5} \epsilon_{\mathrm{F}}  \approx  15 \,\text{MeV}.\end{aligned} \hspace{\stretch{1}}(1.73)

Question: Neutron star ([2], prob 9.5)

Model a neutron star as an ideal Fermi gas of neutrons at T = 0 moving in the gravitational field of a heavy point mass M at the center. Show that the pressure P obeys the equation

\begin{aligned}\frac{dP}{dr} = - \gamma M \frac{\rho(r)}{r^2},\end{aligned} \hspace{\stretch{1}}(1.74)

where \gamma is the gravitational constant, r is the distance from the center, and \rho(r) is the density which only depends on distance from the center.


In the grand canonical scheme the pressure for a Fermion system is given by

\begin{aligned}\beta P V &= \ln Z_{\mathrm{G}} \\ &= \ln \prod_\epsilon \sum_{n = 0}^1 \left( z e^{-\beta \epsilon}  \right)^n \\ &= \sum_\epsilon \ln \left( 1 + z e^{-\beta \epsilon}  \right).\end{aligned} \hspace{\stretch{1}}(1.75)

The kinetic energy of the particle is adjusted by the gravitational potential

\begin{aligned}\epsilon &= \epsilon_\mathbf{k}- \frac{\gamma m M}{r} \\ &= \frac{\hbar^2 \mathbf{k}^2}{2m}- \frac{\gamma m M}{r}.\end{aligned} \hspace{\stretch{1}}(1.76)

Differentiating eq. 1.75 with respect to the radius, we have

\begin{aligned}\beta V \frac{\partial {P}}{\partial {r}} &= -\beta \frac{\partial {\epsilon}}{\partial {r}}\sum_\epsilon \frac{z e^{-\beta \epsilon}}{ 1 + z e^{-\beta \epsilon} } \\ &= -\beta\left( \frac{\gamma m M}{r^2}  \right)\sum_\epsilon \frac{1}{ z^{-1} e^{\beta \epsilon} + 1} \\ &= -\beta\left( \frac{\gamma m M}{r^2}  \right)\left\langle{{N}}\right\rangle.\end{aligned} \hspace{\stretch{1}}(1.77)

Noting that \left\langle{{N}}\right\rangle m/V is the average density of the particles, presumed radial, we have

\begin{aligned}\boxed{\frac{\partial P}{\partial r} &= -\frac{\gamma M}{r^2} \frac{m \left\langle N \right\rangle}{V} \\ &= -\frac{\gamma M}{r^2} \rho(r).}\end{aligned} \hspace{\stretch{1}}(1.0.78)


[1] M. Abramowitz and I.A. Stegun. \emph{Handbook of mathematical functions with formulas, graphs, and mathematical tables}, volume 55. Dover publications, 1964.

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

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

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , , , , , , , , , , | Leave a Comment »

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

Posted by peeterjoot on March 3, 2013

In A compilation of notes, so far, for ‘PHY452H1S Basic Statistical Mechanics’ I posted a link this compilation of statistical mechanics course notes.

That compilation now all of the following too (no further updates will be made to any of these) :

February 28, 2013 Rotation of diatomic molecules

February 28, 2013 Helmholtz free energy

February 26, 2013 Statistical and thermodynamic connection

February 24, 2013 Ideal gas

February 16, 2013 One dimensional well problem from Pathria chapter II

February 15, 2013 1D pendulum problem in phase space

February 14, 2013 Continuing review of thermodynamics

February 13, 2013 Lightning review of thermodynamics

February 11, 2013 Cartesian to spherical change of variables in 3d phase space

February 10, 2013 n SHO particle phase space volume

February 10, 2013 Change of variables in 2d phase space

February 10, 2013 Some problems from Kittel chapter 3

February 07, 2013 Midterm review, thermodynamics

February 06, 2013 Limit of unfair coin distribution, the hard way

February 05, 2013 Ideal gas and SHO phase space volume calculations

February 03, 2013 One dimensional random walk

February 02, 2013 1D SHO phase space

February 02, 2013 Application of the central limit theorem to a product of random vars

January 31, 2013 Liouville’s theorem questions on density and current

January 30, 2013 State counting

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , | 1 Comment »

PHY452H1S Basic Statistical Mechanics. Problem Set 3: State counting

Posted by peeterjoot on February 10, 2013

[Click here for a PDF of this post with nicer formatting]


This is an ungraded set of answers to the problems posed.

Question: Surface area of a sphere in d-dimensions

Consider a d-dimensional sphere with radius R. Derive the volume of the sphere V_d(R) and its surface area S_d(R) using S_d(R) = dV_d(R)/dR.


Let’s start with some terminology and notation, borrowing from [2]

  1. 1-sphere : 2D circle, with “volume” V_2
  2. 2-sphere : 3D sphere, with volume V_3
  3. 3-sphere : 4D Euclidean hypersphere, with “volume” V_4
  4. 4-sphere : 5D Euclidean hypersphere, with “volume” V_5
  5. \cdots

To calculate the volume, we require a parameterization allowing for expression of the volume element in an easy to integrate fashion. For the 1-sphere, we can use the usual circular and spherical coordinate volume elements

\begin{aligned}V_2(R) = 4 \int_0^R r dr \int_0^{\pi/2} d\theta = 4 \frac{R^2}{2} \frac{\pi}{2} = \pi R^2\end{aligned} \hspace{\stretch{1}}(1.0.1a)

\begin{aligned}V_3(R) = 8 \int_0^R r^2 dr \int_0^{\pi/2} \sin\theta d\theta \int_0^{\pi/2} d\phi= 8 \frac{R^3}{3} (1) \frac{\pi}{2}= \frac{4}{3} \pi R^3\end{aligned} \hspace{\stretch{1}}(1.0.1b)

Here, to simplify the integration ranges, we’ve calculated the “volume” integral for just one quadrant or octet of the circle and sphere respectively, as in (Fig 1)

Fig1: Integrating over just one quadrant of circle or octet of sphere

How will we generalize this to higher dimensions? To calculate the volume elements in a systematic fashion, we can introduce a parameterization and use a Jacobian to change from Cartesian coordinates. For the 1-volume and 2-volume cases those parameterizations were the familiar

\begin{aligned}\begin{aligned}x_1 &= r \cos\theta \\ x_0 &= r \sin\theta\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.2a)

\begin{aligned}\begin{aligned}x_2 &= r \cos\theta \\ x_1 &= r \sin\theta \cos\phi \\ x_0 &= r \sin\theta \sin\phi\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.2b)

Some reflection shows that this generalizes nicely. Let’s use shorthand

\begin{aligned}C_k = \cos\theta_k\end{aligned} \hspace{\stretch{1}}(1.0.3a)

\begin{aligned}S_k = \sin\theta_k,\end{aligned} \hspace{\stretch{1}}(1.0.3b)

and pick V_5, say, a dimension bigger than the 2D or 3D cases that we can do by inspection, we can parameterize with

\begin{aligned}\begin{aligned}x_4 &= r C_4 \\ x_3 &= r S_4 C_3 \\ x_2 &= r S_4 S_3 C_2 \\ x_1 &= r S_4 S_3 S_2 C_1 \\ x_0 &= r S_4 S_3 S_2 S_1.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.4)

Our volume element

\begin{aligned}dV_5 = \underbrace{\frac{\partial(x_4, x_3, x_2, x_1, x_0)}{\partial(r, \theta_4, \theta_3, \theta_2, \theta_1)} }_{\equiv J_5}dr d\theta_4 d\theta_3 d\theta_2 d\theta_1\end{aligned} \hspace{\stretch{1}}(1.0.5)

That Jacobian is

\begin{aligned}J_5 = \left\lvert\begin{array}{lllll}C_4 & - r S_4 & \phantom{-} 0 & \phantom{-} 0 & \phantom{-} 0 \\ S_4 C_3 & \phantom{-} r C_4 C_3 & - r S_4 S_3 & \phantom{-} 0 & \phantom{-} 0 \\ S_4 S_3 C_2 & \phantom{-} r C_4 S_3 C_2 & \phantom{-} r S_4 C_3 C_2 & - r S_4 S_3 S_2 & \phantom{-} 0 \\ S_4 S_3 S_2 C_1 & \phantom{-} r C_4 S_3 S_2 C_1 & \phantom{-} r S_4 C_3 S_2 C_1 & \phantom{-} r S_4 S_3 C_2 C_1 & - r S_4 S_3 S_2 S_1 \\ S_4 S_3 S_2 S_1 & \phantom{-} r C_4 S_3 S_2 S_1 & \phantom{-} r S_4 C_3 S_2 S_1 & \phantom{-} r S_4 S_3 C_2 S_1 & \phantom{-} r S_4 S_3 S_2 C_1 \end{array}\right\rvert=r^4 S_4^3 S_3^2 S_2^1 S_1^0\left\lvert\begin{array}{lllll} \begin{array}{l} C_4 \end{array} \begin{array}{l} \phantom{S_3 S_2 S_1} \end{array} & \begin{array}{l} -S_4 \end{array} \begin{array}{l} \phantom{S_3 S_2 S_1} \end{array} & \begin{array}{l} \phantom{-}0 \end{array} & \begin{array}{l} \phantom{-}0 \end{array} & \begin{array}{l} \phantom{-}0 \end{array}\\ \begin{array}{l}S_4 \\ S_4 \\ S_4 \\ S_4 \end{array}\boxed{\begin{array}{l} C_3 \\ S_3 C_2 \\ S_3 S_2 C_1 \\ S_3 S_2 S_1 \end{array} }&\begin{array}{l}\phantom{-}C_4 \\ \phantom{-}C_4 \\ \phantom{-}C_4 \\ \phantom{-}C_4 \end{array}\boxed{\begin{array}{l}C_3 \\ S_3 C_2 \\ S_3 S_2 C_1 \\ S_3 S_2 S_1 \end{array}}&\begin{array}{l} - S_3 \\ \phantom{-} C_3 C_2 \\ \phantom{-} C_3 S_2 C_1 \\ \phantom{-} C_3 S_2 S_1 \end{array}&\begin{array}{l}\phantom{-} 0 \\ - S_2 \\ \phantom{-} C_2 C_1 \\ \phantom{-} C_2 S_1 \end{array}&\begin{array}{l}\phantom{-}0 \\ \phantom{-}0 \\ - S_1 \\ \phantom{-}C_1 \end{array}\end{array}\right\rvert\end{aligned} \hspace{\stretch{1}}(1.0.6)

Observe above that the cofactors of both the 1,1 and the 1,2 elements, when expanded along the first row, have a common factor. This allows us to work recursively

\begin{aligned}J_5 = r^4 S_4^3 S_3^2 S_2^1 S_1^0(C_4^2 - - S_4^2)\left\lvert\begin{array}{llll}C_3 & - S_3 & \phantom{-} 0 & \phantom{-} 0 \\ S_3 C_2 & \phantom{-} C_3 C_2 & - S_2 & \phantom{-} 0 \\ S_3 S_2 C_1 & \phantom{-} C_3 S_2 C_1 & \phantom{-} C_2 C_1 & - S_1 \\ S_3 S_2 S_1 & \phantom{-} C_3 S_2 S_1 & \phantom{-} C_2 S_1 & \phantom{-} C_1 \end{array}\right\rvert=r S_4^3 J_4\end{aligned} \hspace{\stretch{1}}(1.0.7)

Similarly for the 4D volume

\begin{aligned}J_4 = r^3 S_3^2 S_2^1 S_1^0\left\lvert\begin{array}{llll}C_3 & - S_3 & \phantom{-} 0 & \phantom{-} 0 \\ S_3 C_2 & \phantom{-} C_3 C_2 & - S_2 & \phantom{-} 0 \\ S_3 S_2 C_1 & \phantom{-} C_3 S_2 C_1 & \phantom{-} C_2 C_1 & - S_1 \\ S_3 S_2 S_1 & \phantom{-} C_3 S_2 S_1 & \phantom{-} C_2 S_1 & \phantom{-} C_1 \end{array}\right\rvert= r^3 S_3^2 S_2^1 S_1^0 (C_2^2 + S_2^2)\left\lvert\begin{array}{lll}\phantom{-} C_2 & - S_2 & \phantom{-} 0 \\ \phantom{-} S_2 C_1 & \phantom{-} C_2 C_1 & - S_1 \\ \phantom{-} S_2 S_1 & \phantom{-} C_2 S_1 & \phantom{-} C_1 \end{array}\right\rvert= r S_3^2 J_3\end{aligned} \hspace{\stretch{1}}(1.0.8)

and for the 3D volume

\begin{aligned}J_3 = r^2 S_2^1 S_1^0\left\lvert\begin{array}{lll}\phantom{-} C_2 & - S_2 & \phantom{-} 0 \\ \phantom{-} S_2 C_1 & \phantom{-} C_2 C_1 & - S_1 \\ \phantom{-} S_2 S_1 & \phantom{-} C_2 S_1 & \phantom{-} C_1 \end{array}\right\rvert= r^2 S_2^1 S_1^0 (C_2^2 + S_2^2)\left\lvert\begin{array}{ll}C_1 & - S_1 \\ S_1 & \phantom{-} C_1 \end{array}\right\rvert= r S_2^1 J_2,\end{aligned} \hspace{\stretch{1}}(1.0.9)

and finally for the 2D volume

\begin{aligned}J_2= r S_1^0\left\lvert\begin{array}{ll}C_1 & - S_1 \\ S_1 & \phantom{-} C_1 \end{array}\right\rvert= r S_1^0.\end{aligned} \hspace{\stretch{1}}(1.0.10)

Putting all the bits together, the “volume” element in the n-D space is

\begin{aligned}dV_n = dr d\theta_{n-1} \cdots d\theta_{1} r^{n-1} \prod_{k=0}^{n-2} (\sin\theta_{k+1})^{k}, \end{aligned} \hspace{\stretch{1}}(1.0.11)

and the total volume is

\begin{aligned}V_n(R) = 2^n \int_0^R r^{n-1} dr \prod_{k=0}^{n-2} \int_0^{\pi/2} d\theta \sin^k \theta= 2^n \frac{R^{n}}{n}\prod_{k=0}^{n-2} \int_0^{\pi/2} d\theta \sin^k \theta= \frac{4 R^2}{n} (n-2) 2^{n-2} \frac{R^{n-2}}{n-2}\prod_{k=n-3}^{n-2} \int_0^{\pi/2} d\theta \sin^k \theta\prod_{k=0}^{n-4} \int_0^{\pi/2} d\theta \sin^k \theta=4 R^2 \frac{n-2}{n} V_{n-2}(R)\int_0^{\pi/2} d\theta \sin^{n-2} \theta\int_0^{\pi/2} d\theta \sin^{n-3} \theta=4 R^2 \frac{n-2}{n} V_{n-2}(R)\frac{\pi}{2 (n-2)}.\end{aligned} \hspace{\stretch{1}}(1.0.12)

Note that a single of these sine power integrals has a messy result (see: statMechProblemSet3.nb)

\begin{aligned}\int_0^{\pi/2} d\theta \sin^{k} \theta=\frac{\sqrt{\pi } \Gamma \left(\frac{k+1}{2}\right)}{2 \Gamma \left(\frac{k}{2}+1\right)},\end{aligned} \hspace{\stretch{1}}(1.0.13)

but the product is simple, and that result, computed in statMechProblemSet3.nb, is inserted above, providing an expression for the n-D volume element as a recurrence relation

\begin{aligned}\boxed{V_{n} = \frac{2 \pi}{n} R^2 V_{n-2}(R)}\end{aligned} \hspace{\stretch{1}}(1.0.14)

With this recurrence relation we can find the volume V_{2n} or V_{2n+1} in terms of V_2 and V_3 respectively. For the even powers we have

\begin{aligned}\begin{array}{lll}V_{2(2)} &= \frac{2 \pi R^2}{4} V_{2} &= \frac{2^1 \pi R^2}{2(2)} V_{2} \\ V_{2(3)} &= \frac{(2 \pi R^2)^2}{(6)(4)} V_{2} &= \frac{2^2 (\pi R^2)^2}{2^2 (3)(2)} V_{2} \\ V_{2(4)} &= \frac{(2 \pi R^2)^3}{(8)(6)(4)} V_{2} &= \frac{2^3 (\pi R^2)^3}{2^3 (4)(3)(2)} V_{2} \end{array}\end{aligned} \hspace{\stretch{1}}(1.0.15)

\begin{aligned}V_{2(n)} = \frac{2^n (\pi R^2)^{n-1}}{(2n)!!} V_{2} = \frac{(\pi R^2)^{n-1}}{(n)!} V_{2} = \frac{(\pi R^2)^{n-1}}{(n)!} \pi R^2\end{aligned} \hspace{\stretch{1}}(1.0.16)

This gives us a closed form expression for the even powers for n \ge 2

\begin{aligned}\boxed{V_{2(n)} = \frac{(\pi R^2)^{n}}{n!}.}\end{aligned} \hspace{\stretch{1}}(1.0.17)

Observe that this also conveniently gives V_2 = \pi R^2, so is actually valid for n \ge 1. Now for the odd powers

\begin{aligned}\begin{array}{lll}V_{2(2)+1} &= \frac{2 \pi R^2}{5} V_{3} &= 3 \frac{2^1 \pi R^2}{5!!} V_{3} \\ V_{2(3)+1} &= 3 \frac{(2 \pi R^2)^2}{7 5 3} V_{3} &= 3 \frac{2^2 (\pi R^2)^2}{7!!} V_{2} \\ V_{2(4)+1} &= 3 \frac{(2 \pi R^2)^3}{(9)(7)(5)(3)} V_{2} &= 3 \frac{2^3 (\pi R^2)^3}{9!!} V_{2} \end{array}\end{aligned} \hspace{\stretch{1}}(1.0.18)

\begin{aligned}V_{2(n)+1} = 3 \frac{2^{n-1} (\pi R^2)^{n-1}}{(2n +1)!!} V_{2} = 3 \frac{2^{n-1} (\pi R^2)^{n-1}}{(2n +1)!!} \frac{4}{3} \pi R^3\end{aligned} \hspace{\stretch{1}}(1.0.19)

So, for n \ge 2 we have

\begin{aligned}\boxed{V_{2(n)+1} = \frac{2^{n+1} \pi^n R^{2 n + 1}}{(2n +1)!!}.}\end{aligned} \hspace{\stretch{1}}(1.0.20)

As with the even powered expression 1.0.17 we see that this is also good for n = 1, yielding V_3 = 4 \pi R^2/3 as required for the 3D spherical volume.

The even and odd power expressions don’t look quite different on the surface, but can be put into a consistent form, when expressed in terms of the gamma function.

For the even powers, using a substitution 2 n = m we have for even values of m \ge 2

\begin{aligned}V_{m} = \frac{\pi^{m/2} R^m}{\Gamma(n + 1)}= \frac{\pi^{m/2} R^m}{\Gamma(m/2 + 1)}.\end{aligned} \hspace{\stretch{1}}(1.0.21)

For the even powers, with the help of [1] we find

\begin{aligned}(2 n + 1)!! = \frac{2^{n+1} \Gamma\left(n + \frac{3}{2}\right)}{\sqrt{\pi}}.\end{aligned} \hspace{\stretch{1}}(1.0.22)

This gives us

\begin{aligned}V_{2(n)+1} = \frac{\pi^{n + 1/2} R^{2 n + 1}}{\Gamma\left(n + \frac{3}{2}\right)}= \frac{\pi^{n + 1/2} R^{2 n + 1}}{\Gamma\left(\frac{(2 n + 1) + 2}{2}\right)}.\end{aligned} \hspace{\stretch{1}}(1.0.23)

Writing m = 2 n + 1, or n = (m - 1)/2 we have for odd values of m \ge 3 a match with the even powered expression of 1.0.21

\begin{aligned}\boxed{V_{m} = \frac{ \pi^{m/2} R^{m} }{ \Gamma\left( m/2 + 1 \right)}.}\end{aligned} \hspace{\stretch{1}}(1.0.24)

We’ve shown that this is valid for any dimension m \ge 2.

Tabulating some values of these for n \in [2, 7] we have respectively

\begin{aligned}\pi r^2,\frac{4 \pi r^3}{3},\frac{\pi ^2 r^4}{2},\frac{8 \pi ^2 r^5}{15},\frac{\pi ^3 r^6}{6},\frac{16 \pi ^3 r^7}{105},\end{aligned} \hspace{\stretch{1}}(1.0.25)

The only task left is computation of the surface area. That comes by inspection and is

\begin{aligned}S_{m}(R) =\frac{d V_m(R)}{d R}= \frac{\pi^{m/2} m R^{m-1}}{\Gamma\left( m/2 + 1 \right)}.\end{aligned} \hspace{\stretch{1}}(1.0.26)

Again for m \in [2, 7] we have

\begin{aligned}2 \pi r,4 \pi r^2,2 \pi ^2 r^3,\frac{8 \pi ^2 r^4}{3},\pi ^3 r^5,\frac{16 \pi ^3 r^6}{15}.\end{aligned} \hspace{\stretch{1}}(1.0.27)

Question: State counting – polymer

A typical protein is a long chain molecule made of very many elementary units called amino acids – it is an example of a class of such macromolecules called polymers. Consider a protein made N amino acids, and assume each amino acid is like a sphere of radius a. As a toy model assume that the protein configuration is like a random walk with each amino acid being one “step”, i.e., the center-to-center vector from one amino acid to the next is a random vector of length 2a and ignore any issues with overlapping spheres (so-called “excluded volume” constraints). Estimate the spatial extent of this protein in space. Typical proteins assume a compact form in order to be functional. In this case, taking the constraint of nonoverlapping spheres, estimate the radius of such a compact protein assuming it has an overall spherical shape with fully packed amino acids (ignore holes in the packing, and use only volume ratios to make this estimate). With N = 300 and a \approx 5 \text{\AA}, estimate these sizes for the random walk case as well as the compact globular case.


We are considering a geometry like that of (Fig 2), depicted in two dimensions for ease of illustration.

Fig2: Touching “spheres”

From the geometry, if \mathbf{c}_k is the vector to the center of the kth sphere, we have for some random unit vector \hat{\mathbf{r}}_k

\begin{aligned}\mathbf{c}_{k+1} = \mathbf{c}_k + 2 a \hat{\mathbf{r}}_k\end{aligned} \hspace{\stretch{1}}(1.0.28)

Proceeding recursively, writing d_N, and n = N -1 > 0, we have for the difference of the positions of the first and last centers of the chain

\begin{aligned}d_N = \left\lvert {\mathbf{c}_N - \mathbf{c}_1} \right\rvert = 2 a \left\lvert { \sum_{k = 1}^{n} \hat{\mathbf{r}}_k} \right\rvert= 2 a \left( \sum_{k,m = 1}^{n} \hat{\mathbf{r}}_k \cdot \hat{\mathbf{r}}_m \right)^{1/2}= 2 a \left( n + 2 \sum_{1 \le k < m \le n} \hat{\mathbf{r}}_k \cdot \hat{\mathbf{r}}_m \right)^{1/2}= 2 a \sqrt{n} \left( 1 + \frac{2}{n} \sum_{1 \le k < m \le n} \hat{\mathbf{r}}_k \cdot \hat{\mathbf{r}}_m \right)^{1/2}\end{aligned} \hspace{\stretch{1}}(1.0.29)

The \hat{\mathbf{r}}_k‘s clearly cannot be completely random since we have a constraint that \hat{\mathbf{r}}_k \cdot \hat{\mathbf{r}}_{k+1} > -\cos\pi/3, or else two adjacent spheres will overlap. There will also be overlapping constraints for longer portions of the chain that are harder to express. We are ignoring both such constraints, and seek the ensemble average of all systems of the form 1.0.29.

Employing random azimuthal and polar angular variables \theta_k, and \phi_k, we have

\begin{aligned}\hat{\mathbf{r}}_k = \begin{bmatrix}\sin\theta_k \cos\phi_k \\ \sin\theta_k \sin\phi_k \\ \cos\theta_k,\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.30)

so that the average polymer length is

\begin{aligned}\begin{aligned}\left\langle{{d_N }}\right\rangle &=2 a \sqrt{n} \left( \frac{1}{{(2 \pi)(\pi)}} \right)^{n} \int_{\theta_j \in [0, \pi]}d\theta_1 d\theta_2 \cdots d\theta_n \int_{\phi_j \in [0, 2 \pi]}d\phi_1 d\phi_2 \cdots d\phi_n \times \\ &\quad \left( 1 + \frac{2}{n} \sum_{1 \le k < m \le n} \sin\theta_k \cos\phi_k \sin\theta_m \cos\phi_m + \sin\theta_k \sin\phi_k \sin\theta_m \sin\phi_m + \cos\theta_k \cos\theta_m \right)^{1/2}\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.31)

Observing that even \int \sqrt{1 + a \cos\theta} d\theta is an elliptic integral, we don’t have any hope of evaluating this in closed form. However, to first order, we have

\begin{aligned}\begin{aligned}d_N&\approx2 a \sqrt{n} \left( \frac{1}{{(2 \pi)(\pi)}} \right)^{n} \int_{\theta_j \in [0, \pi]}d\theta_1 d\theta_2 \cdots d\theta_n \int_{\phi_j \in [0, 2 \pi]}d\phi_1 d\phi_2 \cdots d\phi_n \times \\ &\quad \left( 1 + \frac{1}{n} \sum_{1 \le k < m \le n} \sin\theta_k \cos\phi_k \sin\theta_m \cos\phi_m + \sin\theta_k \sin\phi_k \sin\theta_m \sin\phi_m + \cos\theta_k \cos\theta_m \right) \\ &=2 a \sqrt{n} \left( \frac{1}{{(2 \pi)(\pi)}} \right)^{n} \left( 2 \pi^2 \right)^n \\ &\quad +2 a \sqrt{n} \left( \frac{1}{{(2 \pi)(\pi)}} \right)^{n} \left( 2 \pi^2 \right)^{n - 2} \left( \frac{1}{{2}} n(n+1) \right) \frac{1}{{n}} \times \\ &\quad \int_0^{\pi}d\theta \int_0^{\pi}d\theta'\int_0^{2 \pi}d\phi\int_0^{2 \pi}d\phi'\left( \sin\theta \cos\phi \sin\theta' \cos\phi' + \sin\theta \sin\phi \sin\theta' \sin\phi' + \cos\theta \cos\theta' \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.32)

The \int_0^{2 \pi} \cos\phi integrals kill off the first term, the \int_0^{2 \pi} \sin\phi integrals kill of the second, and the \int_0^\pi \cos\theta integral kill of the last term, and we are left with just

\begin{aligned}d_N \approx 2 a \sqrt{N-1}.\end{aligned} \hspace{\stretch{1}}(1.0.33)

Ignoring the extra 2 \times a of the end points, and assuming that for large N we have \sqrt{N-1} \approx \sqrt{N}, the spatial extent of the polymer chain is

\begin{aligned}\boxed{d_N \approx 2 a \sqrt{N}.}\end{aligned} \hspace{\stretch{1}}(1.0.34)

Spherical packing

Assuming the densest possible spherical packing, a face centered cubic [3], as in (Fig 3), we see that the density of such a spherical packing is

Fig3: Element of a face centered cubic

\begin{aligned}\eta_{\mathrm{FCC}} =\frac{\left( 8 \times \frac{1}{{8}} + 6 \times \frac{1}{{2}} \right) \frac{4}{3} \pi r^3}{ \left( \sqrt{8 r^2} \right)^3 }= \frac{\pi}{\sqrt{18}}.\end{aligned} \hspace{\stretch{1}}(1.0.35)

With a globular radius of R and an atomic radius of a, and density \eta we have

\begin{aligned}\eta \frac{4}{3} \pi R^3 = N \frac{4}{3} \pi a^3,\end{aligned} \hspace{\stretch{1}}(1.0.36)

so that the globular radius R is

\begin{aligned}\boxed{R_N(\eta) = a \sqrt[3]{\frac{N}{\eta}}.}\end{aligned} \hspace{\stretch{1}}(1.0.37)

Some numbers

With N = 300 and a \approx 5 \text{\AA}, and ignoring spaces (i.e. \eta = 1, for a non-physical infinite packing), our globular diameter is approximately

\begin{aligned}2 \times 5 \text{\AA} \sqrt[3]{300} \approx 67 \text{\AA}.\end{aligned} \hspace{\stretch{1}}(1.0.38)

This is actually not much different than the maximum spherical packing of an FCC lattice, which results a slightly larger globular cluster diameter

\begin{aligned}2 \times 5 \text{\AA} \sqrt[3]{300 \sqrt{18}/\pi} \approx 74 \text{\AA}.\end{aligned} \hspace{\stretch{1}}(1.0.39)

Both however, are much less than the end to end length of the random walk polymer chain

\begin{aligned}2 (5 \text{\AA}) \sqrt{300} \approx 173 \text{\AA}.\end{aligned} \hspace{\stretch{1}}(1.0.40)

Question: State counting – spins

Consider a toy model of a magnet where the net magnetization arises from electronic spins on each atom which can point in one of only two possible directions – Up/North or Down/South. If we have a system with N spins, and if the magnetization can only take on values \pm 1 (Up = +1, Down = -1), how many configurations are there which have a total magnetization m, where m = (N_\uparrow - N_\downarrow)/N (note that N_\uparrow + N_\downarrow = N)? Simplify this result assuming N \gg 1 and a generic m (assume we are not interested in the extreme case of a fully magnetized system where m = \pm 1). Find the value of the magnetization m for which the number of such microscopic states is a maximum. For N = 20, make a numerical plot of the number of states as a function of the magnetization (note: -1 \le m \le 1) without making the N \gg 1 assumption.


For the first couple values of N, lets enumerate the spin sample spaces, their magnetization.

N = 1

  1. \uparrow : m = 1
  2. \downarrow, m = -1

N = 2

  1. \uparrow \uparrow : m = 1
  2. \uparrow \downarrow : m = 0
  3. \downarrow \uparrow, m = 0
  4. \downarrow \downarrow, m = -1

N = 3

  1. \uparrow \uparrow \uparrow : m = 1
  2. \uparrow \uparrow \downarrow : m = 1/3
  3. \uparrow \downarrow \uparrow : m = 1/3
  4. \uparrow \downarrow \downarrow : m = -1/3
  5. \downarrow \uparrow \uparrow : m = 1/3
  6. \downarrow \uparrow \downarrow : m = -1/3
  7. \downarrow \downarrow \uparrow : m = -1/3
  8. \downarrow \downarrow \downarrow : m = 1

The respective multiplicities for N = 1,2,3 are \{1\}, \{1, 2, 1\}, \{1, 3, 3, 1\}. It’s clear that these are just the binomial coefficients. Let’s write for the multiplicities

\begin{aligned}g(N, m) = \binom{N}{i(m)}\end{aligned} \hspace{\stretch{1}}(1.0.41)

where i(m) is a function that maps from the magnetization values m to the integers [0, N]. Assuming

\begin{aligned}i(m) = a m + b,\end{aligned} \hspace{\stretch{1}}(1.0.42)

where i(-1) = 0 and i(1) = N, we solve

\begin{aligned}\begin{aligned}a (-1) + b &= 0 \\ a (1) + b &= N,\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.43)


\begin{aligned}i(m) = \frac{N}{2}(m + 1)\end{aligned} \hspace{\stretch{1}}(1.0.44)


\begin{aligned}g(N, m) = \binom{N}{\frac{N}{2}(1 + m)} = \frac{N!}{\left(\frac{N}{2}(1 + m)\right)!\left(\frac{N}{2}(1 - m)\right)!}\end{aligned} \hspace{\stretch{1}}(1.0.45)


\begin{aligned}\begin{aligned}2 m &= N_{\uparrow} - N_{\downarrow} \\ N &= N_{\uparrow} + N_{\downarrow},\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.46)

we see that this can also be written

\begin{aligned}\boxed{g(N, m) = \frac{N!}{(N_{\uparrow})!(N_{\downarrow})!}}\end{aligned} \hspace{\stretch{1}}(1.0.47)

Simplification for large N

Using Stirlings approximation

\begin{aligned}\begin{aligned}\ln g(N, m) &= \ln N! - \ln N_{\downarrow}! - \ln N_{\uparrow}! \\ &\approx \left( N + \frac{1}{{2}} \right) \ln N - \not{{N}} + \not{{ \frac{1}{{2}} \ln 2 \pi }} \\ &\quad-\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln N_{\uparrow} + \not{{N_{\uparrow}}} - \not{{\frac{1}{{2}} \ln 2 \pi}} \\ &\quad-\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln N_{\downarrow} + \not{{N_{\downarrow}}} - \frac{1}{{2}} \ln 2 \pi \\ &=\left( N_{\uparrow} + N_{\downarrow} + \frac{1}{{2}} + \underbrace{\frac{1}{{2}}}_{\text{add}} \right) \ln N -\frac{1}{{2}} \ln 2 \pi - \underbrace{\frac{1}{{2}} \ln N}_{\text{then subtract}} \\ &\quad -\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln N_{\uparrow} \quad -\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln N_{\downarrow} \\ &=-\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln N_{\uparrow}/N-\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln N_{\downarrow}/N- \frac{1}{{2}} \ln 2 \pi N \\ &=-\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln \frac{1}{{2}}( 1 + m )-\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln \frac{1}{{2}} ( 1 - m )- \frac{1}{{2}} \ln 2 \pi N \\ &\approx\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln 2+\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln 2 - \frac{1}{{2}} \ln 2 \pi N \\ &\quad -\left( N_{\uparrow} + \frac{1}{{2}} \right) \left(m - \frac{1}{{2}} m^2 \right) \\ &\quad -\left( N_{\downarrow} + \frac{1}{{2}} \right) \left( -m - \frac{1}{{2}} m^2 \right) \\ &= (N + 1) \ln 2 - \frac{1}{{2}} \ln 2 \pi N+ m ( N_{\downarrow} - N_{\uparrow} )+ \frac{1}{{2}} m^2 ( N + 1 ) \\ &= (N + 1) \ln 2 - \frac{1}{{2}} \ln 2 \pi N- N m^2+ \frac{1}{{2}} m^2 ( N + 1 ) \\ &\approx(N + 1) \ln 2 - \frac{1}{{2}} \ln 2 \pi N- \frac{1}{{2}} m^2 N\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.48)

This gives us for large N

\begin{aligned}g(N, m) \approx 2^N \sqrt{\frac{2}{ \pi N}} e^{ -m^2 N/2 }\end{aligned} \hspace{\stretch{1}}(1.0.49)

For N = 20 this approximation and the exact expression are plotted in (Fig 4).

Fig4: Distribution of number of configurations for N = 20 magnets as a function of magnetization

With the large scales of this extremely peaked function, no visible difference can be seen. That difference does exist, and is plotted in (Fig 5)

Fig5: N = 20 differences between the exact binomial expression and the Gaussian approximation


[1] M. Abramowitz and I.A. Stegun. \emph{Handbook of mathematical functions with formulas, graphs, and mathematical tables}, volume 55. Dover publications, 1964.

[2] Wikipedia. N-sphere — Wikipedia, The Free Encyclopedia, 2013\natexlab{a}. URL [Online; accessed 26-January-2013].

[3] Wikipedia. Sphere packing — wikipedia, the free encyclopedia, 2013\natexlab{b}. URL [Online; accessed 31-January-2013].

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , | Leave a Comment »

PHY452H1S Basic Statistical Mechanics. Lecture 6: Volumes in phase space. Taught by Prof. Arun Paramekanti

Posted by peeterjoot on January 29, 2013

[Click here for a PDF of this post with nicer formatting]


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

Liouville’s theorem

We’ve looked at the continuity equation of phase space density

\begin{aligned}0 = \frac{\partial {\rho}}{\partial {t}} + \sum_{i_\alpha} \left(\frac{\partial {}}{\partial {p_{i_\alpha}}} \left( \rho \dot{p}_{i_\alpha} \right) + \frac{\partial {\left( \rho \dot{x}_{i_\alpha} \right) }}{\partial {x_{i_\alpha}}}\right)\end{aligned} \hspace{\stretch{1}}(1.2.1)

which with

\begin{aligned}\frac{\partial {\dot{p}_{i_\alpha}}}{\partial {p_{i_\alpha}}} + \frac{\partial {\dot{x}_{i_\alpha}}}{\partial {x_{i_\alpha}}} = 0\end{aligned} \hspace{\stretch{1}}(1.2.2)

led us to Liouville’s theorem

\begin{aligned}\\ boxed{\frac{d{{\rho}}}{dt}(x, p, t) = 0}.\end{aligned} \hspace{\stretch{1}}(1.2.3)

We define Ergodic, meaning that with time, as you wait for t \rightarrow \infty, all available phase space will be covered. Not all systems are necessarily ergodic, but the hope is that all sufficiently complicated systems will be so.

We hope that

\begin{aligned}\rho(x, p, t \rightarrow \infty) \implies \frac{\partial {\rho}}{\partial {t}} = 0 \qquad \mbox{in steady state}\end{aligned} \hspace{\stretch{1}}(1.2.4)

In particular for \rho = \text{constant}, we see that our continuity equation 1.2.1 results in 1.2.2.

For example in a SHO system with a cyclic phase space, as in (Fig 1).

Fig 1: Phase space volume trajectory


\begin{aligned}\left\langle{{A}}\right\rangle = \frac{1}{{\tau}} \int_0^\tau dt A( x_0(t), p_0(t) ),\end{aligned} \hspace{\stretch{1}}(1.2.5)

or equivalently with an ensemble average, imagining that we are averaging over a number of different systems

\begin{aligned}\left\langle{{A}}\right\rangle = \frac{1}{{\tau}} \int dx dp A( x, p ) \underbrace{\rho(x, p)}_{\text{constant}}\end{aligned} \hspace{\stretch{1}}(1.2.6)

If we say that

\begin{aligned}\rho(x, p) = \text{constant} = \frac{1}{{\Omega}},\end{aligned} \hspace{\stretch{1}}(1.2.7)

so that

\begin{aligned}\left\langle{{A}}\right\rangle = \frac{1}{{\Omega}} \int dx dp A( x, p ) \end{aligned} \hspace{\stretch{1}}(1.2.8)

then what is this constant. We fix this by the constraint

\begin{aligned}\int dx dp \rho(x, p) = 1\end{aligned} \hspace{\stretch{1}}(1.2.9)

So, \Omega is the allowed “volume” of phase space, the number of states that the system can take that is consistent with conservation of energy.

What’s the probability for a given configuration. We’ll have to enumerate all the possible configurations. For a coin toss example, we can also ask how many configurations exist where the sum of “coin tosses” are fixed.

A worked example: Ideal gas calculation of \Omega

  • N gas atoms at phase space points \mathbf{x}_i, \mathbf{p}_i
  • constrained to volume V
  • Energy fixed at E.

\begin{aligned}\Omega(N, V, E) = \int_V d\mathbf{x}_1 d\mathbf{x}_2 \cdots d\mathbf{x}_N \int d\mathbf{p}_1 d\mathbf{p}_2 \cdots d\mathbf{p}_N \delta \left(E - \frac{\mathbf{p}_1^2}{2m}- \frac{\mathbf{p}_2^2}{2m}\cdots- \frac{\mathbf{p}_N^2}{2m}\right)=\underbrace{V^N}_{\text{Real space volume, not N dimensional ``volume''}} \int d\mathbf{p}_1 d\mathbf{p}_2 \cdots d\mathbf{p}_N \delta \left(E - \frac{\mathbf{p}_1^2}{2m}- \frac{\mathbf{p}_2^2}{2m}\cdots- \frac{\mathbf{p}_N^2}{2m}\right)\end{aligned} \hspace{\stretch{1}}(1.10)

With \gamma defined implicitly by

\begin{aligned}\frac{d\gamma}{dE} = \Omega\end{aligned} \hspace{\stretch{1}}(1.3.11)

so that with Heavyside theta as in (Fig 2).

\begin{aligned}\Theta(x) = \left\{\begin{array}{l l}1 & \quad x \ge 0 \\ 0 & \quad x < 0\end{array}\right.\end{aligned} \hspace{\stretch{1}}(1.0.12a)

\begin{aligned}\frac{d\Theta}{dx} = \delta(x),\end{aligned} \hspace{\stretch{1}}(1.0.12b)

Fig 2: Heavyside theta


we have

\begin{aligned}\gamma(N, V, E) = V^N \int d\mathbf{p}_1 d\mathbf{p}_2 \cdots d\mathbf{p}_N \Theta \left(E - \sum_i \frac{\mathbf{p}_i^2}{2m}\right)\end{aligned} \hspace{\stretch{1}}(1.0.13)

In three dimensions (p_x, p_y, p_z), the dimension of momentum part of the phase space is 3. In general the dimension of the space is 3N. Here

\begin{aligned}\int d\mathbf{p}_1 d\mathbf{p}_2 \cdots d\mathbf{p}_N \Theta \left(E - \sum_i \frac{\mathbf{p}_i^2}{2m}\right),\end{aligned} \hspace{\stretch{1}}(1.0.14)

is the volume of a “sphere” in 3N– dimensions, which we found in the problem set to be

\begin{aligned}V_{m} = \frac{ \pi^{m/2} R^{m} }{ \Gamma\left( m/2 + 1 \right)}.\end{aligned} \hspace{\stretch{1}}(1.0.15a)

\begin{aligned}\Gamma(x) = \int_0^\infty dy e^{-y} y^{x-1}\end{aligned} \hspace{\stretch{1}}(1.0.15b)

\begin{aligned}\Gamma(x + 1) = x \Gamma(x) = x!\end{aligned} \hspace{\stretch{1}}(1.0.15c)

Since we have

\begin{aligned}\mathbf{p}_1^2 + \cdots \mathbf{p}_N^2 \le 2 m E\end{aligned} \hspace{\stretch{1}}(1.0.16)

the radius is

\begin{aligned}\text{radius} = \sqrt{ 2 m E}.\end{aligned} \hspace{\stretch{1}}(1.0.17)

This gives

\begin{aligned}\gamma(N, V, E) = V^N \frac{ \pi^{3 N/2} ( 2 m E)^{3 N/2}}{\Gamma( 3N/2 + 1) }= V^N \frac{2}{3N} \frac{ \pi^{3 N/2} ( 2 m E)^{3 N/2}}{\Gamma( 3N/2 ) },\end{aligned} \hspace{\stretch{1}}(1.0.17)


\begin{aligned}\Omega(N, V, E) = V^N \pi^{3 N/2} ( 2 m E)^{3 N/2 - 1} \frac{2 m}{\Gamma( 3N/2 ) }\end{aligned} \hspace{\stretch{1}}(1.0.19)

This result is almost correct, and we have to correct in 2 ways. We have to fix the counting since we need an assumption that all the particles are indistinguishable.

  • Indistinguishability. We must divide by N!.
  • \Omega is not dimensionless. We need to divide by h^{3N}, where h is Plank’s constant.

In the real world we have to consider this as a quantum mechanical system. Imagine a two dimensional phase space. The allowed points are illustrated in (Fig 3).

Fig 3: Phase space volume adjustment for the uncertainty principle


Since \Delta x \Delta p \sim \hbar, the question of how many boxes there are, we calculate the total volume, and then divide by the volume of each box. This sort of handwaving wouldn’t be required if we did a proper quantum mechanical treatment.

The corrected result is

\begin{aligned}\boxed{\Omega_{\mathrm{correct}} = \frac{V^N}{N!} \frac{1}{{h^{3N}}} \frac{( 2 \pi m E)^{3 N/2 }}{E} \frac{1}{\Gamma( 3N/2 ) }}\end{aligned} \hspace{\stretch{1}}(1.0.20)

To come

We’ll look at entropy

\begin{aligned}\underbrace{S}_{\text{Entropy}} = \underbrace{k_{\mathrm{B}}}_{\text{Boltzmann's constant}} \ln \underbrace{\Omega_{\mathrm{correct}}}_{\text{phase space volume (number of configurations)}}\end{aligned} \hspace{\stretch{1}}(1.0.21)

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , | Leave a Comment »

Hypersphere volume calculation the easy way

Posted by peeterjoot on January 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.)]


In problem set three I solved the hypersphere volume problem the hard way. I did it the way that I thought was obvious, starting with the spherical coordinate hypervolume element for an Euclidean space. That build on the volume element calculation I’d previously done for four dimensional Euclidean and Hyperbolic spaces in [1]. This time I avoided any use of Geometric Algebra and wrote out the volume element directly using a Jacobian transformation instead of a wedge product (that leads to the same Jacobian). I then proceeded to integrate the volume element, come up with a recurrence relation for the volume, solve that for even and odd powers, and find a common expression using Gamma functions that was correct for even and odd powers. It was a labourious process, but satisfying since I’d previously tried this calculation and not succeeded.

As and after I did this calculation, which I’ll post later, I figured there had to have been an easier way. In the midterm prep reading of section 5.5 of [2] I found just that method. It’s done there in about six lines, using a trick that seemed really sneaky! I was left with a frustrated feeling, wondering how on earth somebody figured out to do it that way.

After a bit of reflection, I see that the trick is a very sensible approach. I’ll outline that here for 2 and 3 dimensions to show the technique, and the reader can generalize if desired.

The trick

Switching to spherical or circular coordinates when there is radial symmetry isn’t anything that we could describe as trickery. For example, with r^2 = x^2 + y^2 in a two dimensional problem or r^2 = x^2 + y^2 + z^2 in a three dimensional problem, it’s clear that we’d do the following respectively if we were evaluating an integral

\begin{aligned}\iint f(r) dx dy = \iint f(r) r dr d\theta = 2 \pi \int f(r) r dr\end{aligned} \hspace{\stretch{1}}(1.0.5a)

\begin{aligned}\iiint f(r) dx dy = \iiint f(r) r^2 dr \sin d\theta d\theta d\phi = 4 \pi \int f(r) r^2 dr\end{aligned} \hspace{\stretch{1}}(1.0.5b)

In fact, for f(r) = 1 these give us the area and volume of the circle and sphere respectively.

So, what’s the trick? The first part is the observation that the “area” of a “volume” for the circle and the sphere are found by the derivatives of the “volumes”

\begin{aligned}\frac{d}{dr} \pi r^2 = 2 \pi r\end{aligned} \hspace{\stretch{1}}(1.0.2a)

\begin{aligned}\frac{d}{dr} \frac{4 \pi}{3} r^3 = 4 \pi r^2.\end{aligned} \hspace{\stretch{1}}(1.0.2b)

I recall being suprised that this is the case back in high school calculus. When I lost marks for not having the formula for the surface area of a sphere memorized on some test, my calculus teacher told me this title tidbit.

Back then this wasn’t obvious to me, and I complained that it was true only for the perfect symmetry of a sphere. For example, the derivative of an volume of a cube doens’t give the surface area of a cube (d x^3/dx = 3 x^2 \ne 6 x^2).

Once we believe or assume that the surface “area” of a hypervolume is the derivative of the volume, we can proceed with the trick. That trick is to express the volume in terms of an unknown constant. For example, for the circle and sphere the generalized “volume”s are respectively

\begin{aligned}V_2 = B r^2\end{aligned} \hspace{\stretch{1}}(1.0.3a)

\begin{aligned}V_3 = C r^3\end{aligned} \hspace{\stretch{1}}(1.0.3b)

The perimeter of the circle and surface area of the sphere are then

\begin{aligned}A_2 = 2 B r\end{aligned} \hspace{\stretch{1}}(1.0.4a)

\begin{aligned}A_3 = 3 C r^2\end{aligned} \hspace{\stretch{1}}(1.0.4b)

So, if we want to calculate integrals of the form we can write

\begin{aligned}\iint f(r) dx dy = 2 B \int f(r) r dr\end{aligned} \hspace{\stretch{1}}(1.0.5a)

\begin{aligned}\iiint f(r) dx dy = 3 C \int f(r) r^2 dr.\end{aligned} \hspace{\stretch{1}}(1.0.5b)

The essence of the trick is to do such an integral in both Cartesian coordinates to get the left hand sides, and then do the radial right hand side integrals. Comparing these provides the constants B and C and thus completes the “volume” formulas for the circle and sphere.

The function chosen for this integral in the text was a Gaussian exponential f(r) = e^{-r^2/2}, something that is strictly radial, and can be integrated over all space. For the 2D case, we’ll integrate

\begin{aligned}\iint e^{-(x^2 +y^2)/2} dx dy = (\sqrt{2 \pi})^2= 2 B \int_0^\infty dr r e^{-r^2/2}= -2 B {\left.{e^{-r^2/2}}\right\vert}_{0}^{\infty} = 2 B.\end{aligned} \hspace{\stretch{1}}(1.0.5b)

We find that B = \pi, so the 2D spherical “volume” (the area of the circle) is V = \pi r^2.

For the 3D sphere, we have

\begin{aligned}\iiint e^{-(x^2 +y^2 + z^2)/2} dx dy = (\sqrt{2 \pi})^3= 3 C \int_0^\infty dr r^2 e^{-r^2/2}= 3 C \int_0^\infty dr e^{-r^2/2}= 3 C \sqrt{2 \pi}.\end{aligned} \hspace{\stretch{1}}(1.0.5b)

So we have for the volume of the 3D sphere, V = 4 \pi r^3/3 as expected. The same idea can be extended to higher dimensional spheres. The text does the even values of N. Treating both even and odd values, I’d expect this to yield the result I calculated with the Jacobian volume element method

\begin{aligned}V_{m} = \frac{ \pi^{m/2} R^{m} }{   \Gamma\left( m/2 + 1 \right)}.\end{aligned} \hspace{\stretch{1}}(1.0.5b)


[1] Peeter Joot. Exploring physics with Geometric Algebra, chapter {Spherical and hyperspherical parametrization}. URL

[2] S.K. Ma. Statistical Mechanics. World Scientific, 1985. ISBN 9789971966072. URL

Posted in Math and Physics Learning. | Tagged: , , , , , , | 3 Comments »