\documentclass[11pt,fleqn]{article}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% template that does not use Revtex4
%%% but allows special fonts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Please use this template.
%%% Edit it using e.g. Notepad
%%% Ignore the header (do not change it)
%%% Process the file in the Latex site 
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Page setup
\topmargin -1.5cm      
\oddsidemargin -0.04cm   
\evensidemargin -0.04cm  
\textwidth 16.59cm
\textheight 24cm 
\setlength{\parindent}{0cm} 
\setlength{\parskip}{0cm} 
\usepackage{graphicx}

% Fonts
\usepackage{latexsym}
\usepackage{amsmath} 
\usepackage{amssymb} 
\usepackage{bm}



% Math symbols I
\newcommand{\sinc}{\mbox{sinc}}
\newcommand{\const}{\mbox{const}}
\newcommand{\trc}{\mbox{trace}}
\newcommand{\intt}{\int\!\!\!\!\int }
\newcommand{\ointt}{\int\!\!\!\!\int\!\!\!\!\!\circ\ }
\newcommand{\ar}{\mathsf r}
\newcommand{\im}{\mbox{Im}}
\newcommand{\re}{\mbox{Re}}

% Math symbols II
\newcommand{\eexp}{\mbox{e}^}
\newcommand{\bra}{\left\langle}
\newcommand{\ket}{\right\rangle}

% Mass symbol
\newcommand{\mass}{\mathsf{m}} 
\newcommand{\Mass}{\mathsf{M}} 

% More math commands
\newcommand{\tbox}[1]{\mbox{\tiny #1}}
\newcommand{\bmsf}[1]{\bm{\mathsf{#1}}} 
\newcommand{\amatrix}[1]{\begin{matrix} #1 \end{matrix}} 
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}

% Other commands
\newcommand{\hide}[1]{}
\newcommand{\drawline}{\begin{picture}(500,1)\line(1,0){500}\end{picture}}
\newcommand{\bitem}{$\bullet$ \ \ \ }
\newcommand{\Cn}[1]{\begin{center} #1 \end{center}}
\newcommand{\mpg}[2][1.0\hsize]{\begin{minipage}[b]{#1}{#2}\end{minipage}}
\newcommand{\Dn}{\vspace*{3mm}}

% Figures
\newcommand{\putgraph}[2][0.30\hsize]{\includegraphics[width=#1]{#2}}

% heading
\newcommand{\exnumber}[1]{\newcommand{\exnum}{#1}}
\newcommand{\heading}[1]{\begin{center} {\Large {\bf Ex\exnum:} #1} \end{center}}
\newcommand{\auname}[1]{\begin{center} {\bf Submitted by:} #1 \end{center}}


\begin{document} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\exnumber{5024}
\heading{Pressure of Lenard Jones gas} 
\auname{Rami Rahimi}


{\bf The problem:}

\Dn

A gas of $N$ particles is confined in a box 
of volume $V$ at temprature of $T$.
The two-body interaction between the particles is 
given by the Lenard Jones expression:

\[ u(r)\ \ =\ \ \frac {a} {r^{12}} - \frac{b}{r^6} \]

Note that this interaction is characterized by 
a length scale $r_0$ and an energy scale $\epsilon_0$ 
that correspond to the position and the depth of the potential.    


\Dn

(a) Find an expression for the pressure via the Virial theorem, assuming that the moments $\langle r^n \rangle_T$ are known.

(b) Using the Virial expansion, find an explicit expression for the pressure assuming low temperatures.

(c) Using the Virial expansion, find an explicit expression for the pressure assuming high temperatures.

(d) Comparing your answers to items (a) and (c) deduce explicit expressions for the $n=-6$ and for the $n=-12$ moments. Express your result in terms of ${ (V, r_0, \epsilon_0, T) }$. 



\Dn\Dn


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf The solution:}

\Dn

The Lenard Jones potential is:
\begin{align} 
u(r)= 4\epsilon_0\left[\left(\frac{r_0}{r} \right)^{12} - \left(\frac{r_0}{r} \right)^{6} \right]
\end{align}
The minimum of this potential is at $r_m=2^{1/6}r_0$, such that $u(r_m)=-{\epsilon_0}$ and $u''(r_m)={72{\epsilon_0}}/{r_m^2} $, where $\epsilon_0$ and $r_0$ are the energy and length scales respectively.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\Dn

(a) The pressure is given by the Virial theorem:
\begin{align} 
P={\frac{1}{V}}\left[NT - \frac{1}{3}\big \langle {r\cdot\frac{\partial U}{\partial r}}\big \rangle  \right]
\end{align}
The two-body interaction between the particles has the form:
\begin{align} 
U=\sum_{ \langle ij \rangle }4\epsilon_0\left(\left(\frac{r_0}{\|r_i-r_j \|} \right)^{12} - \left(\frac{r_0}{\|r_i-r_j \|} \right)^{6} \right)= \sum_{ \langle ij \rangle }U_{ij}
\end{align}
So, the pressure is:
\begin{align} 
P={\frac{1}{V}}\left[NT - \frac{1}{3}\big \langle {\sum_{ \langle ij \rangle }\left(r_i\cdot\frac{\partial U_{ij}}{\partial r_i}+r_j\cdot\frac{\partial U_{ij}}{\partial r_j}\right) }\big \rangle  \right]  \end{align}
The partial derivatives are done carefully:
\begin{align} 
\frac{\partial U_{ij}}{\partial r_i}\propto\frac{\partial}{\partial r_i} \|r_i-r_j \|^\alpha={\alpha \|r_i-r_j \|^{\alpha-2}}{\left( \vec{r_i}-\vec{r_j} \right)}
\end{align}
Notice that $\frac{\partial U_{ij}}{\partial r_j}=-\frac{\partial U_{ij}}{\partial r_i}$ hence:
\begin{align} 
\left(\vec{r_i}\cdot\frac{\partial U_{ij}}{\partial r_i}+\vec{r_j}\cdot\frac{\partial U_{ij}}{\partial r_j}\right)={\alpha\left(\vec{r_i} -\vec{r_j}\right)\left(\vec{r_i} -\vec{r_j}\right)}{ \|r_i-r_j \|^{\alpha-2}}={\alpha \|r_i-r_j \|^{\alpha}}=\alpha U_{ij}
\end{align}
And finally we get: 
\begin{equation} \label{eq:a}
P={\frac{N}{V}}\Big[T + {4\epsilon_0\,r_0^6N}\left(2r_0^{6} \langle r^{-12} \rangle- \langle r^{-6} \rangle \right) \Big]
\end{equation}
Note that the sum over $ \langle ij \rangle $ has been replaced by the total number of the interactions which is ${N\left(N-1\right)}/{2}  \approx {N^2}/{2}$.

%%%%%%%%%%%%%%%
\Dn

(b)\ The Virial expansion for the pressure is: 
\begin{equation} \label{eq:press}
P=\frac{NT}{V} \bigg[1+a_2\left(\frac{N}{V}\right)+a_3\left(\frac{N}{V}\right)^2+...\bigg]
\end{equation}
The second Virial coefficient is:
\begin{align} 
a_2=-\frac{1}{2}\int {\left(e^{-\beta u\left(r\right)} -1\right)} d^3r=-\frac{1}{2}\int {f\left(r\right)} d^3r 
\end{align}

{\begin{figure}[h!]
  \centering
    \includegraphics[width=0.5\textwidth]{5024a}
  \caption{Comparison between $f(r)$ with different values of $T$ }
 \label{fig:gull}
\end{figure}
}
\Dn

The area where $f(r)$ is negative originates from the repulsive part of $u$ and is bounded by $[0, r_0] $ in the horizontal axis and by $[−1, 0]$ in the vertical axis, while the positive area of $f(r)$ originates from the attraction is unbounded and increases as the temperature $T$ decreases since $f(r_m)\propto \exp{\big(\frac{\epsilon_0}{T}\big)}$ (as shown in Figure 1). Additionally, for low temperatures, the -1 term in the integrand can be ignored and the integral can be calculated by saddle point approximation around the minimum point of $u(r)$. This approximation is valid due the fact that the main contribution to the integral is around $r_m$ since $f(r) \to 0$ for $r\gg r_m$:
\begin{align} 
a_2\approx - \displaystyle\frac{1}{2} \displaystyle\int {e^{-\beta u\left(r_m\right)-\frac{1}{2}\beta u''\left(r_m\right){\left(r-r_m\right)^2}}} d^3r= - \displaystyle\frac{1}{2}e^{-\beta u\left(r_m\right)}4\pi \int{e^{-\frac{1}{2}\beta u''\left(r_m\right){\left(r-r_m\right)^2}}}r^2 dr 
\end{align}
The last integral (up to normalization factor of $\sqrt{2\pi\sigma^2}$) has the form of $\langle r^2 \rangle$ which is the second moment of gaussian with $\sigma^2=(\beta u''(r_m))^{-1}$ and $\langle r \rangle=r_m$, so it can be calculated directly by the variance definition $ \langle r^2 \rangle=\sigma^2+{\langle r \rangle}^2$:
\begin{align} 
\displaystyle a_2=-2\pi e^{-\beta u \left(r_m \right)}r_m^3\sqrt{\frac{\pi T}{36\epsilon_0}}{\left(\frac{T}{72\epsilon_0}+1 \right)}\approx -\frac{r_0^3}{3}\sqrt{\frac{2\pi^3T}{\epsilon_0}}\exp\left({\frac{\epsilon_0}{T}}\right)
\end{align}
The last approximation is valid under the assumption of low $T$. Hence, the pressure is:
\begin{equation} 
P=\frac{NT}{V} \bigg[1 -\frac{r_0^3}{3}\frac{N}{V}\sqrt{\frac{2\pi^3T}{\epsilon_0}}\exp\left({\frac{\epsilon_0}{T}}\right)\bigg]
\end{equation}

%%%%%%%%%%%%%%
\Dn
(c)
In the limit of high temperature i.e. $ \epsilon_0/T\ll1$ the integrand $f\left(r\right)$ can be expanded up to first order in Taylor:
\begin{align} 
a_2=\displaystyle -\frac{1}{2}\int {\left(e^{-\beta u\left(r\right)} -1\right)} d^3r\approx -\frac{1}{2}\int {\Big({-\beta u\left(r\right)}\Big)} d^3r
\end{align}
For the region $r>r_0$ this approximation is valid, but for the limit of small $r$ (i.e. $r\ll r_0$) the approximation deviates from the original function  since $ -\beta u(r) \to -\infty $ while $f(r)\to -1$. Figure 2 demonstrates the difference between $f(r)=\left(e^{-\beta u\left(r\right)} -1\right)$ and $-\beta u(r)$.\\

\begin{figure}[h!]
  \centering
    \includegraphics[width=0.5\textwidth]{5024b}
  \caption{Comparison between $f(r)$ and its Taylor expansion $ -\beta u\left(r\right) $ }
 \label{fig:gull}
\end{figure}





In order to solve this problem, the integral over $f(r)$, meaning $a_2$, is divided into two regions $r>r_0$ and $r<r_0$   such that $a_2=a_{>}+a_{<}$ respectively. \\
For the $r > r_0$ part, the Taylor expansion is valid, and we get:\\
\begin{equation} \label{eq:highr}
 a_{>}=-\frac{1}{2}\int_{r_0}^{\infty}{-\beta u\left(r\right) d^3r}=8\pi\beta\epsilon_0\int_{r_0}^{\infty}{\left[\left(\frac{r_0}{r} \right)^{12} - \left(\frac{r_0}{r} \right)^{6} \right] r^2dr}=-\frac{16\pi}{9}\Big(\frac{\epsilon_0}{T}\Big)r_0^3
\end{equation}
\vspace{3Cm}

In the region $r < r_0$, we have $f(0) = -1$ and $f(r_0) = 0$. For small enough $r$, the function $f(r) \approx -1$, while for $r \approx r_0$, we have $f(r) \approx 0$. Hence we approximate $f(r)$ as a step function with
the width $\tilde{r}$, so that $\beta U(\tilde{r})=1$. For high temperature $\tilde{r}$ is small, such that $\tilde{r}=r_0( {{4\epsilon_0}/{T}})^{{1}/{12}}$. Figure 3 demonstrates this approximation.


{\begin{figure}[h!]
  \centering
    \includegraphics[width=0.5\textwidth]{5024cn}
  \caption{Enlargement of the negative area of $f(r)$}
 \label{fig:gull}
\end{figure}}



Consequently, for this hard-core potential approximation, $a_2$ has the form:
\begin{align} \label{eq:lowr}
a_{<}=\displaystyle -\frac{1}{2}\int { -1 } \ d^3r=\frac{4\pi}{2}\int_{0}^{\tilde{r}}{r^2dr}=\frac{2\pi}{3} \tilde{r}^3=\displaystyle\frac{2\pi}{3}\Big(\frac{4\epsilon_0}{T}\Big)^{{1}/{4}}r_0^3
\end{align}
Placing Eqs~\ref{eq:highr} and~\ref{eq:lowr} into the total $a_2$ in the Virial expansion (Eq~\ref{eq:press}) will give the pressure:
\begin{align} \label{eq:c}
P=\frac{NT}{V} \bigg[1+\frac{2\pi}{3}\Big(\frac{4\epsilon_0}{T}\Big)^{{1}/{4}}r_0^3\left(\frac{N}{V}\right)-\frac{16\pi\epsilon_0}{9T}r_0^3\left(\frac{N}{V}\right)\bigg]
\end{align}


%%%%%%%%%%%%%

(d)
The positive and negative items in Eqs~\ref{eq:a} and~\ref{eq:c} can be compared separately:
\begin{align} 
\displaystyle\frac{8N^2\epsilon_0r_0^{12}}{V} \langle r^{-12} \rangle=\frac{2\pi}{3}\Big(\frac{4\epsilon_0}{T}\Big)^{{1}/{4}}r_0^3\left(\frac{N^2T}{V^2}\right)  \Longrightarrow \langle r^{-12} \rangle=\frac{\pi}{3Vr_0^9}\Big(\frac{T}{4\epsilon_0}\Big)^{{3}/{4}}
\end{align}
\begin{align} 
\displaystyle\frac{4N^2\epsilon_0r_0^{6}}{V} \langle r^{-6} \rangle=\frac{16\pi\epsilon_0}{9T}r_0^3\left(\frac{N^2T}{V^2}\right) \Longrightarrow \langle r^{-6} \rangle=\frac{4\pi}{9Vr_0^3} 
\end{align}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document} 
