Homework #2

Homework #2#

1. Solid or Gas Earth?#

In class, we defined the quantity:

\[\Gamma = \frac{E_\mathrm{Coul}}{kT}\]

as a measure of the importance of the Coulomb energy to the thermal energy. If \(\Gamma \gg 1\), then we are a solid. We got the separation between nuclei from the number density and the temperature from the Virial theorem and got:

\[\Gamma = \frac{Z^2 e^2}{A^{4/3} m_u^{4/3} G M^{2/3}}\]

Now for the earth, we can take \(A\) and \(Z\) to just be for iron

A = 56
Z = 26

For the constants, we’ll use a library, unyt to get them along with the units

import unyt

e = unyt.qe.in_cgs()
m_u = (1.0 * unyt.amu).in_cgs()
G = unyt.G.in_cgs()
M = unyt.Mearth.in_cgs()

Notice that these all have units

m_u

unyt_quantity(1.66053892e-24, 'g')

unyt_quantity(-4.80320467e-10, 'statC')

gamma = (Z * e)**2 / ((A * m_u)**(4./3.) * G * M**(2./3.))
gamma

unyt_quantity(167.13713567, '(dimensionless)')

So this tells us that gamma is dimensionless!

\[\Gamma = 167\]

Now, really that we need to be careful with Coulomb’s law and units. In CGS, the energy is:

\[E_\mathrm{Coul} = \frac{q_1 q_2}{r}\]

but in SI, it is:

\[E_\mathrm{Coul} = \frac{1}{4\pi \epsilon_0} \frac{q_1 q_2}{r}\]

Let’s redo that calculation but using SI units:

e = unyt.qe
m_u = 1 *unyt.amu
G = unyt.G
M = unyt.Mearth

gamma = (Z * e)**2 / ((A * m_u)**(4./3.) * G * M**(2./3.))
gamma

unyt_quantity(2.20215509e-17, 'C**2*s**2/(amu**(1/3)*kg**(2/3)*m**3)')

Notice that the units don’t work now! Let’s add in the \(1/(4\pi \epsilon_0)\)

import math

eps0 = unyt.eps_0

gamma = (Z * e)**2 / ((A * m_u)**(4./3.) * G * M**(2./3.)) / (4 * math.pi * eps0)
gamma

unyt_quantity(167.13713567, '(dimensionless)')

2. Ideal gas#

from sympy import init_session
init_session(use_latex="mathjax")

IPython console for SymPy 1.13.3 (Python 3.11.10-64-bit) (ground types: python)

These commands were executed:
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at https://docs.sympy.org/1.13.3/

ni, m, k, T = symbols("n_I m k_B T", real=True, positive=True)
p, v = symbols("p v", real=True)
n_I = symbols("n_I", real=True)
g = symbols("g", real=True, positive=True)
h = symbols("h", real=True, positive=True)
E_0 = symbols("E_0", real=True)
mu = symbols("mu", real=True)

a. Finding Maxwell#

We start with

n = 4 * pi * p**2 * g / h**3 * exp(mu/(k*T)) * exp(-E_0/(k*T)) * exp(-p**2/(2*m*k*T))
n

\[\displaystyle \frac{4 \pi g p^{2} e^{- \frac{E_{0}}{T k_{B}}} e^{\frac{\mu}{T k_{B}}} e^{- \frac{p^{2}}{2 T k_{B} m}}}{h^{3}}\]

now let’s integrate over p

n_I_int = integrate(n, (p, 0, oo))
n_I_int

\[\displaystyle \frac{2 \sqrt{2} \pi^{\frac{3}{2}} T^{\frac{3}{2}} g k_{B}^{\frac{3}{2}} m^{\frac{3}{2}} e^{- \frac{E_{0}}{T k_{B}}} e^{\frac{\mu}{T k_{B}}}}{h^{3}}\]

This is equal to \(n_I\). We want to substitute this back into \(n(p)\). Let’s solve for \(g\) in terms of \(n_I\) and then use that as a way to substitute this back.

g_solve = solve(Eq(n_I_int, n_I), g)[0]

g_solve

\[\displaystyle \frac{\sqrt{2} h^{3} n_{I} e^{\frac{E_{0} - \mu}{T k_{B}}}}{4 \pi^{\frac{3}{2}} T^{\frac{3}{2}} k_{B}^{\frac{3}{2}} m^{\frac{3}{2}}}\]

n = simplify(n.subs(g, g_solve, strict=False))
n

\[\displaystyle \frac{\sqrt{2} n_{I} p^{2} e^{- \frac{p^{2}}{2 T k_{B} m}}}{\sqrt{\pi} T^{\frac{3}{2}} k_{B}^{\frac{3}{2}} m^{\frac{3}{2}}}\]

Now we see that we have the normal Maxwell-Boltzmann distribution.

b. Ideal gas pressure#

The pressure integral is simply:

\[P = \frac{1}{3} \int_0^\infty n(p) v p dp\]

with \(v = p / m\).

P = integrate(Rational(1, 3) * n * p**2 / m, (p, 0, oo))
P

\[\displaystyle T k_{B} n_{I}\]

We see that we get the familiar ideal gas law!

c. Energy#

The energy is:

\[\rho e = \int_0^\infty n(p) E(p) dp\]

with \(E(p) = p^2 / (2m)\)

rhoe = integrate(n * p**2 / (2 *m), (p, 0, oo))
rhoe

\[\displaystyle \frac{3 T k_{B} n_{I}}{2}\]

So we see that this is \(3/2 P\).