Stirling's approximation and paramagnetism

Stirling's approximation

How large is \(n!\), and how accurate is Stirling's approximation? To refresh our memory about Matlab, and to get a feel for just how HUGE some of these numbers are, let's make a plot of \(\log(n!)\).

n = 100;
r = 10:n;
lognfac = log(factorial(r));
xlabel('n'); ylabel('log(n!)');

How many digits are there in 100! ?

ans =


158 digits! By contrast, \(10^{100}\) only has 100 digits. So 100! is much, much bigger than a googol. The total number of atoms in the whole universe is estimated to be \(10^{80}\) or there about. Thus 100! is roughly the square of the number of particles in universe!

Stirling's approximation is the following:

\(\log(n!) = n \log(n) - n + \frac{1}{2} \log(2\pi n) + O(1/n)\,.\)

Often only the first two terms are good enough for applications. Occasionally the third term is also relevant. How accurate is this? Let's plot the residual error

a1 = r.*log(r);
a2 = a1 - r;
a3 = a2 + .5*log(2*pi*r);

approx1 = abs( lognfac - a1 );
approx2 = abs( lognfac - a2 );
approx3 = abs( lognfac - a3 );

title('Absolute error for k terms in Stirling approximation');
xlabel('n'); ylabel('|log(n!) - a_n|');

Looks like the first approximation is not so good. What is the percent error for each of these approximations?

relapprox1 = approx1./lognfac;
relapprox2 = approx2./lognfac;
relapprox3 = approx3./lognfac;

title('Relative error for k terms in Stirling approximation');
xlabel('n'); ylabel('|1 - a_k / log(n!)|');

As you can see, the percent error is already quite small for the two-term approximation. By the third, it's totally negligable.

Flipping coins

Suppose we have a coin that lands heads with p=.5 (a fair coin), and we flip it 50 times.

  1. How many possible outcomes are there?
  2. How many ways are there of getting 25 heads and 25 tails?
  3. What is the probability of getting exactly 25/25, in any order?
  4. What is the probability of getting 30/20, in any order? 40/10?
  5. Plot a graph of Pr(k) where k is the number of heads.
  6. Repeat that plot for several values of p.
nchoosek(n,k) % number of ways to get k heads in n flips
ans*p^k*(1-p)^(n-k) % probability of k head in n flips

% Define a function for the probability of k heads
binom = @(n,k,p) nchoosek(n,k)*p^k*(1-p)^(n-k);

% Here is a combined plot for two values, p=0.1 and p=0.5.
for p=[.1,.5]
for k = 0:50
    b(k+1) = binom(n,k,p);
hold on
title('Probability of heads in 50 weighted coin flips (p=0.1, 0.5)');
axis([0 50 0 0.25]);
xlabel('k'); ylabel('Pr(k)');
hold off
ans =


ans =


Two-level Paramagnet

A paramagnet is a system which only has magnetization in the presence of an external magnetic field. This is in contrast with a ferromagnet which has strong interactions between the atoms (or molecules, etc.) that can cause it to retain a magnetization even in the presence of zero external field. A simple model for a paramagnet is a collection of \(n\) non-interacting particles with two energy levels each. We are going to consider just a single particle.

The energy levels of a paramagnet are described by the energy function

\(E_{\mu} = -\vec{\mu}\cdot \vec{B} = - \mu B \cos \theta\,.\)

The two-state paramagnet can only be either aligned or antialigned with the magnetic field, so the only allowed energy levels are \(\pm \mu B\). The aligned state has the lower (\(-\)) energy.

  1. What is the partition function for a single particle?

For a single particle, we just sum the exponentials of the energies:

\(Z = \e^{-\mu B/kT}+\e^{+\mu B/kT} = 2 \cosh(\mu B/kT)\,.\)

Now suppose that our collection of paramagnets is at \(T=300K\) (room temperature), so \(kT = .026\) eV. Suppose that \(\mu = 1\) eV/T.

  1. Compute the formula (as a function of the applied field \(\vec{B} = B\, \hat{z}\,\)) for the average energy.
  2. Compute the average magnetic moment
  3. Plot the both of these as a function of \(B\).

We can derive the first formula as follows:

\(\begin{aligned} \avg{E}\, & =(-\mu B) \e^{+\mu B/kT}/Z + (+\mu B) \e^{-\mu B/kT}/Z \\ & = (-\mu B) \frac{\e^{+\mu B/kT} - \e^{-\mu B/kT}}{2 \cosh(\mu B/kT)} \\ & = -\mu B \tanh(\mu B/kT)\,. \end{aligned}\)

The second formula is given by

\(\begin{aligned} \avg{\vec{\mu}\cdot \hat{z}}\, & = (+\mu) \e^{+\mu B/kT}/Z + (-\mu) \e^{-\mu B/kT}/Z \\ & = \mu \tanh(\mu B/kT)\,. \end{aligned}\)

We can plot these as follows:

kT = .026;
mu = 1;
B = -.1:.001:.1;
aveenergy = -mu*B.*tanh(mu*B/kT);
avemag = mu*tanh(mu*B/kT);
title('Average mu and E of 1eV paramagnet vs. B at 300 K');
xlabel('B (Tesla)'); legend('average \mu (eV/T)', 'average E (eV)');

Note, \(.005 T\) is the strength of a typical refrigerator magnet, and \(0.15 T\) is the strength of a sunspot.