In [1]:

```
%pylab inline
```

Populating the interactive namespace from numpy and matplotlib

In [16]:

```
from sympy import *
from pylab import *
```

In [3]:

```
import matplotlib.pyplot as plt
```

In [4]:

```
from sympy.interactive import printing
printing.init_printing(use_latex=True)
from IPython.display import display
a,x,Z = symbols('a,x,Z')
```

The partition function for 1D system is

$Z = \sum_i e^{-\beta E_i}$

The probability on a state $E_i$

$$P(E_i) = \frac{e^{-\beta E_i}}{Z}$$

$$ M = \mu N \left( e^{\beta \mu B} - e^{-\beta \mu B} \right)/Z$$

Write magnetization as a function of $T$ and $T$.

$$M(T,B) = \mu N \frac{ \left( e^{\mu B /(k_B T)} - e^{- \mu B/(k_B T)} \right) } { \left( e^{\mu B /(k_B T)} + e^{- \mu B/(k_B T)} \right) } = \mu N \tanh(\mu B /(k_B T))$$

$$\frac{M(T)}{\bar M} = e^{\bar T/T} - e^{-\bar T/T} $$

in which, $\bar M = \mu N$ .

In [5]:

```
tanh(0)
```

Out[5]:

$$0.0$$

At $T\rightarrow \infty$, we have $\frac{M}{\bar M}$ 1.

In [6]:

```
tanh(inf)
```

Out[6]:

$$1.0$$

Now plot out this result.

In [7]:

```
t = linspace(0, 100, 1000)
mt = tanh(1/t)
fig11 = plt.figure()
axes = fig11.add_axes([0.1, 0.1, 0.8, 0.8]) # left, bottom, width, height (range 0 to 1)
axes.plot(t, mt, 'r')
axes.set_xlabel('T bar')
axes.set_ylabel('M bar')
axes.set_title('M ~ T');
```

-c:2: RuntimeWarning: divide by zero encountered in divide

This plot agrees with our results of the limits.

Define a unit magnetic induction $\bar B = k_B T/\mu$. Write the magnetization as a function of $B$.

$$ \frac{M(B)}{\bar M} = \tanh(B/\bar B) $$

Plot the result.

In [8]:

```
b = linspace(0, 5, 100)
mb = tanh(b)
fig12 = plt.figure()
axes = fig12.add_axes([0.1, 0.1, 0.8, 0.8]) # left, bottom, width, height (range 0 to 1)
axes.plot(b, mb, 'r')
axes.set_xlabel('B/Bbar')
axes.set_ylabel('M/Mbar')
axes.set_title('M ~ B');
```

Energy of the system is

$$ E = - \mu B N \tanh(\mu B /(k_B T)) $$

In [9]:

```
t13e = linspace(0, 5, 100)
e1 = -tanh(t13e)
fig13e = plt.figure()
axes = fig13e.add_axes([0.1, 0.1, 0.8, 0.8]) # left, bottom, width, height (range 0 to 1)
axes.plot(t13e, e1, 'r')
axes.set_xlabel('T/Tbar')
axes.set_ylabel('E/Ebar')
axes.set_title('E ~ T');
```

Heat capacity is

$$ C = \frac{\partial}{\partial T} E = \frac{\mu^2 B^2 N}{k_B} \frac{ (1 - \tanh^2(\mu B/(k_B T)) }{T^2} $$

In [10]:

```
t13 = linspace(0, 5, 100)
c1 = tanh(t13)
fig13 = plt.figure()
axes = fig13.add_axes([0.1, 0.1, 0.8, 0.8]) # left, bottom, width, height (range 0 to 1)
axes.plot(t13, c1, 'r')
axes.set_xlabel('T/Tbar')
axes.set_ylabel('C/Cbar')
axes.set_title('C ~ T');
```

Partition function is

$$ Z = \int _ \theta 2 e^{-\beta E(\theta)} \mathrm d \theta = 2\pi \mathrm {BesselI}(0, \mu B \beta) $$

where $\mathrm{BesselI}$ is the modified Bessel function of the first kind.

Then we can calculate magnetization, which is given by

$$ M = \frac{1}{Z}\int _ \theta 2N\mu \cos\theta e^{-\beta E(\theta)} d \theta $$

The result given by Mathematica is (so bad that I can't make sympy to deal with such integrations.)

$$ M = \mu N \frac{ \mathrm{BesselI}(1, B \beta \mu) }{ \mathrm{BesselI}(0, B \beta \mu) }$$

Energy should be

$$ E = - \mu B N \frac{ \mathrm{BesselI}(1, B \beta \mu) }{ \mathrm{BesselI}(0, B \beta \mu) }$$

Heat capacity is the derivitive of energy over temperature.

$$ C = B^2 \mu^2 N \frac{-2 \mathrm{BesselI}(0, B \beta \mu)^2 + \mathrm{BesselI}(0,\mu\beta B)( \mathrm{BesselI}(0,\mu\beta B) + \mathrm{BesselI}(2,\mu\beta B) ) }{2 k T^2 \mathrm{BesselI}(0, B \beta \mu)^2} $$

Plot out these results using Mathematica.

Energy of a dipole is $E(\theta) = -\mu B \cos\theta$.

Then we can calculate the partition function, which is

$$ Z = \int_\theta \int _ \phi N e^{-\beta E(\theta)}\sin\theta d\theta d \phi = 4\pi N\frac{ \sinh(\mu B \beta)}{\mu B \beta} $$

Magnetization is

$$ M = -\frac{\mu N}{\mu B \beta} + \mu N \coth(\mu B \beta) = \mu N \left( - \frac{\bar B}{B} + \coth \left(\frac{B}{\bar B}\right) \right) $$

$$ \frac{M}{\bar M} = - \frac{\bar B}{B} + \coth \left(\frac{B}{\bar B}\right) $$

Total energy is

$$ E = \frac{1}{Z} \int _ \phi \int _ \theta (-\mu B \cos\theta N) e^{-\mu E(\theta)}\sin\theta d \theta d\phi = - N \left( \frac{1}{\beta} - B \mu \coth(\mu \beta B) \right) $$

Take the derivative with respect to temperature, we have the heat capacity.

$$ C = N \left( 1 - \left(\frac{\bar T}{T} \right)^2 \mathrm{csch}^2 \left( \frac{\bar T}{T} \right) \right) $$

Plot out the results.

$$ \mathrm d F \leq 0 $$

Suppose we have a gas system going through a reversible process, free energy

$$ \mathrm d F = - S\mathrm d T - p \mathrm d V $$

For fixed $T$ and $V$, we have $\mathrm d F = 0$. For irreversible process, it is less than 0.

Chemical potential is the increased Gibbs free energy when one mole particle is added to the system.