# Mixed systems: effect of interaction energy and the solubility gap

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problem 5.58.

So far when we’ve examined the Gibbs free energy and entropy of a mixture of two pure substances, we’ve assumed that the interaction energy is the same between all molecules in the mixture, that is, that a species ${A}$ molecule has the same interaction energy with another ${A}$ molecule as it does with a ${B}$ molecule (and likewise for ${B}$ molecules).

To get an idea of how a difference in the interaction energies affects the phase transition of a mixture, we’ll consider a simple model as follows. We have our usual two species ${A}$ and ${B}$, and the potential interaction energy between two nearest-neighbour ${A}$s is the same as between two ${B}$s, and we’ll call this energy ${u_{0}}$. Between a nearest neighbour ${A}$ and ${B}$, however, the energy is different, say ${u_{AB}}$. Further, suppose each molecule has ${n}$ nearest neighbours, and any given molecule interacts only with its nearest neighbours. For a system composed of unmixed samples of ${A}$ and ${B}$, with total number ${N=N_{A}+N_{B}}$, we can work out the total potential energy as follows.

Any particular molecule has ${n}$ other molecules with which it interacts, so its total interaction energy is ${nu_{0}}$. There are a total of ${N}$ molecules, but if we simply multiply ${nu_{0}}$ by ${N}$, we are counting each interaction twice, so the total interaction energy is

$\displaystyle U_{unmixed}=\frac{1}{2}Nnu_{0} \ \ \ \ \ (1)$

Now suppose we mix the two populations thoroughly. Let ${x}$ be the proportion of molecules that are of species ${B}$, so that ${1-x}$ is the proportion for ${A}$. On average, an ${A}$ molecule now has ${xn}$ nearest neighbours of type ${B}$ and ${\left(1-x\right)n}$ of type ${A}$, so its interaction energy is now

$\displaystyle U_{A}=xnu_{AB}+\left(1-x\right)nu_{0} \ \ \ \ \ (2)$

Similarly for type ${B}$:

$\displaystyle U_{B}=\left(1-x\right)nu_{AB}+xnu_{0} \ \ \ \ \ (3)$

To get the total energy, we multiply 2 by ${N_{A}=\left(1-x\right)N}$ and 3 by ${N_{B}=xN}$, add and divide by 2 to avoid double-counting:

 $\displaystyle U_{mixed}$ $\displaystyle =$ $\displaystyle \frac{1}{2}\left[xnu_{AB}+\left(1-x\right)nu_{0}\right]\left(1-x\right)N+\frac{1}{2}\left[\left(1-x\right)nu_{AB}+xnu_{0}\right]xN\ \ \ \ \ (4)$ $\displaystyle$ $\displaystyle =$ $\displaystyle x\left(1-x\right)nNu_{AB}+\frac{1}{2}\left(\left(1-x\right)^{2}+x^{2}\right)nNu_{0}\ \ \ \ \ (5)$ $\displaystyle$ $\displaystyle =$ $\displaystyle x\left(1-x\right)nNu_{AB}+\frac{1}{2}\left[1-2x\left(1-x\right)\right]nNu_{0} \ \ \ \ \ (6)$

We can now subtract 1 from this to get the energy change resulting from mixing:

 $\displaystyle \Delta U$ $\displaystyle =$ $\displaystyle U_{mixed}-U_{unmixed}\ \ \ \ \ (7)$ $\displaystyle$ $\displaystyle =$ $\displaystyle x\left(1-x\right)nN\left(u_{AB}-u_{0}\right) \ \ \ \ \ (8)$

Depending on the sign of ${u_{AB}-u_{0}}$, a plot of ${\Delta U}$ looks like this:

The red line is for ${u_{AB}-u_{0}>0}$ and the blue line for ${u_{AB}-u_{0}<0}$.

Note that the slopes of ${\Delta U}$ are finite at the endpoints:

$\displaystyle \frac{d\Delta U}{dx}=\left(1-2x\right)nN\left(u_{AB}-u_{0}\right) \ \ \ \ \ (9)$

which is finite when ${x=0}$ or ${x=1}$.

As we saw earlier, the Gibbs free energy of an unmixed system is

$\displaystyle G=\left(1-x\right)G_{A}^{\circ}+xG_{B}^{\circ} \ \ \ \ \ (10)$

and the change in entropy due to mixing is

$\displaystyle \Delta S_{mixing}=-Nk\left[x\ln x+\left(1-x\right)\ln\left(1-x\right)\right] \ \ \ \ \ (11)$

Combining this with 8 we get the Gibbs energy for our system, taking account of both the difference in interaction energy and the entropy of mixing:

$\displaystyle G=\left(1-x\right)G_{A}^{\circ}+xG_{B}^{\circ}+x\left(1-x\right)nN\left(u_{AB}-u_{0}\right)+NkT\left[x\ln x+\left(1-x\right)\ln\left(1-x\right)\right] \ \ \ \ \ (12)$

In order to plot the Gibbs energy, we need to plug in a few numbers. I’ll take ${G_{A}^{\circ}=100\mbox{ J}}$, ${G_{B}^{\circ}=50\mbox{ J}}$, ${N=6.02\times10^{23}}$, ${n=4}$ and ${u_{AB}-u_{0}=10^{-21}\mbox{ J}}$. With these values, we can plot ${G}$ for several temperatures:

From the top down, the red curve is for ${T=50\mbox{ K},}$ dark green for ${T=100\mbox{ K}}$, blue for ${T=144.928\mbox{ K}}$ (the reason for this peculiar value will be explained below), light green for ${T=175\mbox{ K}}$ and yellow for ${T=200\mbox{ K}}$.

Note that the shape of the curve changes from concave down for low temperatures to concave up for high temperatures, with a combination of concavities for the intermediate temperature of ${T=100\mbox{ K}}$ (dark green curve). It’s worth looking at this latter curve in more detail.

Suppose we expand the plot of this curve and draw the straight line (in grey) that is the lowest line that is tangent to the curve at two points:

The grey line represents a system consisting of two unmixed substances, one of which has ${x=0.081}$ (the value of ${x}$ at the left-most tangent point) and the other of which has ${x=0.919}$ (the right-most tangent). For values of ${x}$ between these points, the grey line represents a system which consists of two separate subsystems, one with mostly ${A}$ molecules and the other with mostly ${B}$ molecules. Since the Gibbs energy of the separated system is lower than that of a fully mixed system (represented by the green curve), the system will be found in this relatively unmixed state for these values of ${x}$. This is known as a solubility gap. By the way, note that these two unmixed subsystems are not pure systems, since they both contain molecules of both species. It’s just that one subsystem consists of about 92% ${A}$ molecules while the other consists of about 92% ${B}$ molecules. As we vary ${x}$, the proportion of each of these two subsystems varies, but the composition within each subsystem remains the same: there is no mixing between them.

Returning now to the compound plot of ${G}$ for various temperatures above, we note that the central maximum that is present in the dark green curve gets shallower until at some critical temperature it disappears and then as the temperature is increased further, it turns into a minimum. We can find this critical temperature as follows.

The two tangent points (where the grey line intersects the green curve above) can be found by taking the derivative of 12:

$\displaystyle \frac{dG}{dx}=G_{B}^{\circ}-G_{A}^{\circ}+\left(1-2x\right)nN\left(u_{AB}-u_{0}\right)+NkT\ln\frac{x}{1-x} \ \ \ \ \ (13)$

The tangent points are defined by the condition that their slope is equal to that of an unmixed system, which is ${G_{B}^{\circ}-G_{A}^{\circ}}$, so the tangent points must satisfy

$\displaystyle \left(1-2x\right)nN\left(u_{AB}-u_{0}\right)+NkT\ln\frac{x}{1-x}=0 \ \ \ \ \ (14)$

One solution of this equation is always ${x=0.5}$, since this makes both terms on the LHS zero. This corresponds to the maximum point on the green curve above. To find the two tangent points we’d need to solve the transcendental equation for ${x}$ (which I did using Maple to get the plot above), but we actually don’t need to do that to find the critical temperature. The key is to look at the second derivative at the maximum point. The second derivative is negative for a maximum and positive for a minimum, so the critical temperature (the highest temperature at which a solubility gap exists) is that value of ${T}$ where the second derivative is zero. That is

 $\displaystyle \frac{d^{2}G}{dx^{2}}$ $\displaystyle =$ $\displaystyle -2nN\left(u_{AB}-u_{0}\right)+NkT\left.\left(\frac{1}{x}+\frac{1}{1-x}\right)\right|_{x=0.5}=0\ \ \ \ \ (15)$ $\displaystyle T_{c}$ $\displaystyle =$ $\displaystyle \frac{n\left(u_{AB}-u_{0}\right)}{2k} \ \ \ \ \ (16)$

For the values we used above, ${T_{c}=144.928\mbox{ K}}$, which explains why I chose this value for the blue curve above.

Finally, we can generate a phase diagram which shows the temperature at which a solubility gap appears. From 14, if a tangent occurs at a given value of ${x}$, the corresponding temperature must be

$\displaystyle T=\frac{\left(1-2x\right)n\left(u_{AB}-u_{0}\right)}{k\ln\frac{1-x}{x}} \ \ \ \ \ (17)$

For the parameters we used above, a plot of this is as follows:

For temperatures above the curve, full mixing occurs, while for temperatures below the curve, the system separates into two unmixed subsystems. The maximum on this curve occurs at ${x=0.5}$ and ${T=144.928\mbox{ K}}$, as above.

# Gibbs free energy of a mixture of two ideal gases

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problem 5.56.

To study phase changes of mixtures of substances, rather than pure substances on their own, it’s best to start by looking at the Gibbs free energy of the mixture. If we start with a collection of molecules of two types, ${A}$ and ${B}$, that are initially separated but whose total number ${N_{A}+N_{B}=N}$ is fixed. Since the two populations are separated, the total Gibbs energy is just the sum of the energies of the two individual populations. If population ${B}$ makes up a fraction ${x}$ and ${A}$ a fraction ${1-x}$ of the total, then

$\displaystyle G=\left(1-x\right)G_{A}^{\circ}+xG_{B}^{\circ} \ \ \ \ \ (1)$

What happens if we now mix the two populations, but keep the pressure and temperature constant? The Gibbs energy is defined as

$\displaystyle G\equiv U+PV-TS \ \ \ \ \ (2)$

It’s possible that the energy ${U}$ will change due to interactions between the two species being different than interactions between molecules of the same species. It’s also possible that the volume will change, if the two species either attract or repel each other differently than molecules of the same species. The biggest change, however, is likely to come from a change in entropy, because with the two populations mixed there are now a great many more ways that the molecules can be arranged.

As we saw earlier, if the two substances are ideal gases at the same pressure and temperature and we allow them to mix so that the total volume is unchanged, the change in entropy is

 $\displaystyle \Delta S_{mixing}$ $\displaystyle =$ $\displaystyle -Nk\left[x\ln x+\left(1-x\right)\ln\left(1-x\right)\right]\ \ \ \ \ (3)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -nR\left[x\ln x+\left(1-x\right)\ln\left(1-x\right)\right] \ \ \ \ \ (4)$

where ${n}$ is the total number of moles of both species and ${R}$ is the gas constant.

If we ignore changes in ${U}$ and ${V}$, then the Gibbs energy after mixing is

$\displaystyle G=\left(1-x\right)G_{A}^{\circ}+xG_{B}^{\circ}+nRT\left[x\ln x+\left(1-x\right)\ln\left(1-x\right)\right] \ \ \ \ \ (5)$

The change in entropy looks like this:

The slopes at the two ends are actually infinite, as we can see by taking the derivative of 4:

$\displaystyle \frac{1}{nR}\frac{d\Delta S_{mixing}}{dx}=-\ln x+\ln\left(1-x\right)=\ln\frac{1-x}{x} \ \ \ \ \ (6)$

As ${x\rightarrow0}$, ${\frac{1-x}{x}\rightarrow+\infty}$ so the log also tends to ${+\infty}$. As ${x\rightarrow1}$, ${\frac{1-x}{x}\rightarrow0}$ and the log tends to ${-\infty}$.

A comparison of 1 and 5 looks like this:

Here we’ve taken ${G_{A}^{o}=100}$, ${G_{B}^{o}=75}$, ${n=1}$ and ${T=100}$. The yellow line shows the Gibbs energy from 1, where the two populations are unmixed. The green curve shows 5, where the two populations are mixed. Since the energy for the mixture is lower (because of the increase in entropy), it is the favoured state.

The Gibbs energy also has infinite slopes at ${x=0}$ and 1:

$\displaystyle \frac{dG}{dx}=-G_{A}^{o}+G_{B}^{o}+nRT\ln\frac{x}{1-x} \ \ \ \ \ (7)$

Because the numerator and denominator of the log have swapped places from 6, the slope at ${x=0}$ is now ${-\infty}$ and the slope at ${x=1}$ is now ${+\infty}$.

The infinite slopes show that even if the proportion of one of the species is very small, there is a very large increase in entropy when the two species are allowed to mix, so there is a strong tendency for this to occur.

# Van der Waals equation of state: Maxwell construction

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problems 5.52 – 5.53.

To get a better feel for the van der Waals equation of state, we can look at how the Gibbs free energy varies with the volume and pressure. Schroeder derives an expression for ${G}$ in his equation 5.56, which I will just quote:

$\displaystyle G=-NkT\ln\left(V-Nb\right)+\frac{\left(NkT\right)\left(Nb\right)}{V-Nb}-\frac{2aN^{2}}{V}+c\left(T\right) \ \ \ \ \ (1)$

where ${c\left(T\right)}$ is a function of temperature only, and can be ignored in what follows since we’ll be specifying a particular temperature in each example, so ${c\left(T\right)}$ is effectively a constant in each case.

This equation is easier to work with if we rewrite it in terms of the reduced variables, defined earlier as ${v\equiv V/V_{c}}$, ${t\equiv T/T_{c}}$ and ${p\equiv P/P_{c}}$, with:

 $\displaystyle V_{c}$ $\displaystyle =$ $\displaystyle 3Nb\ \ \ \ \ (2)$ $\displaystyle T_{c}$ $\displaystyle =$ $\displaystyle \frac{8}{27}\frac{a}{bk}\ \ \ \ \ (3)$ $\displaystyle P_{c}$ $\displaystyle =$ $\displaystyle \frac{a}{27b^{2}} \ \ \ \ \ (4)$

If we divide both sides of 1 by ${NkT_{c}}$ and simplify using the above equations, we get

$\displaystyle G_{c}\equiv\frac{G}{NkT_{c}}=-t\ln\left(3v-1\right)+\frac{t}{3v-1}-\frac{9}{4v}+c^{\prime}\left(t\right) \ \ \ \ \ (5)$

where ${c^{\prime}\left(t\right)=c\left(T\right)/NkT_{c}-t\ln\left(Nb\right)}$ is another function that depends only on temperature, so can be ignored in what follows.

To plot ${G_{c}}$ as a function of reduced pressure, we can use the equation derived in the last post:

$\displaystyle p=\frac{8tv^{2}-9v+3}{\left(3v-1\right)v^{2}} \ \ \ \ \ (6)$

We can specify a temperature ${t}$ and then vary ${v}$, using it to calculate ${p}$ and the corresponding ${G_{c}}$ from the last two equations, and thus generate a plot of ${G_{c}}$ versus ${p}$. For ${t=0.95}$ the isotherm looks like this:

If we plot the Gibbs energy, we get:

As we reduce the volume, we move from right to left on the ${pv}$ plot, and the corresponding path on the Gibbs plot is the line coming up from the bottom. At point ${C}$, we can continue to decrease ${v}$, which causes ${G}$ to carry on along the same path above point ${C}$, until it reaches a sharp corner, at which point it decreases until it reaches point ${B}$, which is directly above (that is, it has the same pressure) as point ${C}$. However, note that ${B}$ does not have the minimum Gibbs energy for that pressure, since point ${C}$ is below it. Thus ${B}$ is unstable and does not correspond to a physical state of the substance. Back on the ${pv}$ plot, as we decrease ${v}$ below point ${B}$, the Gibbs path decreases until it reaches the lower sharp corner. This corner corresponds to the minimum between ${A}$ and ${B}$ on the ${pv}$ plot. Again, note that this isn’t the minimum Gibbs energy for that pressure, so this is also unstable. In fact the whole section of the ${pv}$ plot where the pressure decreases with decreasing volume (between the minimum between ${A}$ and ${B}$ and the maximum between ${B}$ and ${C}$) lies on a Gibbs path that is not the minimum energy, so this whole range of volumes is unstable and is not found in the real world.

As we continue to decrease ${v}$ from the minimum point down to point ${A}$, the Gibbs energy returns from the lower sharp corner until it reaches point ${A}$, where the Gibbs energy is equal to the value it had at point ${C}$. Thus points ${A}$ and ${C}$ have the same Gibbs energy at the same pressure, but point ${A}$ has a much lower volume than point ${C}$, so it corresponds to a condensed state (a liquid) while ${C}$ corresponds to a gas.

In fact, this phase transformation occurs for only one pressure, as we can see from the following argument. At constant temperature, the thermodynamic identity states that

$\displaystyle \left(\frac{\partial G}{\partial P}\right)_{T}=V \ \ \ \ \ (7)$

As we go around the triangular loop in the Gibbs plot above, the overall change in ${G}$ is zero, so if we integrate it around the loop, we must get zero:

$\displaystyle 0=\oint dG=\oint\left(\frac{\partial G}{\partial P}\right)_{T}dP=\oint V\;dP \ \ \ \ \ (8)$

If we divide this equation by ${V_{c}P_{c}}$ it changes nothing, since both these quantities depend only on ${T}$, which is held constant. Thus

$\displaystyle \oint v\;dp=0 \ \ \ \ \ (9)$

To do this integral directly requires inverting 6 to get ${v}$ as a function of ${p}$. However, this requires solving a cubic equation in ${v}$, which gives 3 different roots. Just for completeness, I did actually do this using Maple, and if we plot the result we get, for ${t=0.95}$:

The red, blue and green curves represent contributions from each of the three roots. As we go round the triangular loop above, we go from point ${A}$ to ${B}$ to ${C}$. Requiring this integral to be zero is equivalent (as Schroeder shows) to the area bounded by the vertical line and curve between ${B}$ and ${C}$ being equal to the area bounded by the vertical line and curve between ${A}$ and ${B}$. This is known as the Maxwell construction (yes, it’s the same Maxwell that devised Maxwell’s equations in electrodynamics).

The problem is then to find the value ${p_{0}}$ of ${p}$ that makes these two areas equal. Doing this requires moving the yellow line in the plot left or right until the areas are equal. Solving the problem analytically is a bit trickier, of course, since we’re dealing with a cubic equation.

It’s a lot easier to solve if we revert to the traditional ${pv}$ plot above. The area bounded by the vertical line and curve between ${B}$ and ${C}$ is then given by

$\displaystyle A_{BC}=\int_{v_{B}}^{v_{C}}p\;dv-p_{0}\left(v_{C}-v_{B}\right) \ \ \ \ \ (10)$

That is, it’s the area under the ${p}$ curve between ${v_{B}}$ and ${v_{C}}$ minus the area of the rectangle under the horizontal line between the same two volumes. By the same reasoning, the area in the ${vp}$ plot bounded by the vertical line and curve between ${A}$ and ${B}$ is

$\displaystyle A_{AB}=p_{0}\left(v_{B}-v_{A}\right)-\int_{v_{A}}^{v_{B}}p\;dv \ \ \ \ \ (11)$

Surprisingly, the integral turns out to be fairly simple, using Maple:

 $\displaystyle \int p\;dv$ $\displaystyle =$ $\displaystyle \int\frac{8tv^{2}-9v+3}{\left(3v-1\right)v^{2}}dv\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3}{v}+\frac{8}{3}t\ln\left(3v-1\right) \ \ \ \ \ (13)$

The problem now becomes:

1. Specify the temperature.
2. Choose a starting value for the pressure ${p_{0}}$. This can be estimated from the ${pv}$ plot by judging by eye what value of ${p_{0}}$ makes the two areas roughly equal.
3. Solve 6 to find the 3 volumes ${v_{A}}$, ${v_{B}}$ and ${v_{C}}$ for this value of ${p_{0}}$. This requires solving a cubic equation in ${v}$.
4. Calculate the two areas using 10 and 11.
5. Adjust ${p_{0}}$ and repeat from step 3 until we get convergence.

I did this using a Maple procedure, with the best result (for ${t=0.95}$):

$\displaystyle p_{0}=0.8119 \ \ \ \ \ (14)$

As a check, we can compare this with the Gibbs energy plot above and note that the point ${A,C}$ does indeed occur at this pressure.

We can repeat the procedure for ${t=0.8}$ and get the result:

$\displaystyle p_{0}=0.3834 \ \ \ \ \ (15)$

The Gibbs plot for this case is:

Again, we can see that the crossover between the two phases occurs at around ${p_{0}=0.3834}$ so it looks like all is well.

# Surface tension and cloud formation

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problem 5.46.

At the boundary between two phases, there are some molecules that are in a sort of “transition phase”, that is, in neither one phase nor the other. To describe the phase boundary more fully, then, we need to add a term to the Gibbs free energy for these boundary molecules. It turns out that the thickness of the boundary layer is more or less constant, independent of the area of the boundary. The extra Gibbs free energy is therefore proportional to the area of the boundary:

$\displaystyle G_{boundary}=\sigma A \ \ \ \ \ (1)$

where the proportionality coefficient ${\sigma}$ is called the surface tension. It represents the minimum work (per unit area) required to reshape a portion of a substance into a form with a larger surface area. For liquid water at ${20^{\circ}\mbox{ C}}$, ${\sigma=0.073\mbox{ J m}^{-2}}$.

As an example, suppose we have a spherical droplet of water of radius ${r}$ containing ${N_{l}}$ molecules surrounded by ${N-N_{l}}$ molecules of water vapour. If we write the Gibbs free energy in terms of the chemical potential, the Gibbs energy of this system, initially neglecting surface tension, is

$\displaystyle G=N_{l}\mu_{l}+\left(N-N_{l}\right)\mu_{g} \ \ \ \ \ (2)$

where the subscripts ${l}$ and ${g}$ stand for ‘liquid’ and ‘gas’.

Given the volume per molecule ${v_{l}}$ in the liquid, we can write

 $\displaystyle v_{l}$ $\displaystyle =$ $\displaystyle \frac{4\pi r^{3}}{3N_{l}}\ \ \ \ \ (3)$ $\displaystyle G$ $\displaystyle =$ $\displaystyle -\frac{4\pi r^{3}}{3v_{l}}\left(\mu_{g}-\mu_{l}\right)+N\mu_{g} \ \ \ \ \ (4)$

Adding in the energy due to surface tension gives

$\displaystyle G=-\frac{4\pi r^{3}}{3v_{l}}\left(\mu_{g}-\mu_{l}\right)+N\mu_{g}+4\pi r^{2}\sigma \ \ \ \ \ (5)$

Equilibrium values of the sphere’s radius ${r}$ are obtained when ${\partial G/\partial r=0}$, so that there is no change in Gibbs energy when ${r}$ is changed by an amount ${dr}$. This gives

 $\displaystyle \frac{\partial G}{\partial r}$ $\displaystyle =$ $\displaystyle -\frac{4\pi r^{2}}{v_{l}}\left(\mu_{g}-\mu_{l}\right)+8\pi r\sigma=0\ \ \ \ \ (6)$ $\displaystyle r_{c}$ $\displaystyle =$ $\displaystyle \frac{2\sigma v_{l}}{\mu_{g}-\mu_{l}} \ \ \ \ \ (7)$

In order for ${r_{c}>0}$, we need ${\mu_{g}>\mu_{l}}$. The Gibbs energy for this case looks qualitatively like this [units are abitrary; it’s just the shape of the graph we’re interested in]:

Because ${G}$ is a maximum at the equilibrium point, this is an unstable equilibrium, meaning that if ${r}$ is nudged slightly lower, the tendency is for the droplet to continue to decrease in size due to evaporation, while if ${r}$ is nudged slightly higher, the tendency is for the droplet to continue to increase in size.

If ${\mu_{g}<\mu_{l}}$ we get:

There is a minimum at ${r=0}$, but this isn’t a physical solution.

To convert 7 into a form involving the relative humidity, we can start with the thermodynamic identity for the Gibbs energy:

$\displaystyle dG=-S\;dT+V\;dP+\mu\;dN \ \ \ \ \ (8)$

At constant temperature and number we have

$\displaystyle dG=N\;d\mu=V\;dP \ \ \ \ \ (9)$

So for the liquid and gas phases, a slight change in the vapour pressure ${P}$ results in

 $\displaystyle d\mu_{g}$ $\displaystyle =$ $\displaystyle \frac{V_{g}}{N_{g}}dP\ \ \ \ \ (10)$ $\displaystyle$ $\displaystyle =$ $\displaystyle v_{g}dP\ \ \ \ \ (11)$ $\displaystyle d\mu_{l}$ $\displaystyle =$ $\displaystyle v_{l}dP \ \ \ \ \ (12)$

Therefore the change in the difference between the chemical potentials is given by

 $\displaystyle d\left(\mu_{g}-\mu_{l}\right)$ $\displaystyle =$ $\displaystyle \left(v_{g}-v_{l}\right)dP\ \ \ \ \ (13)$ $\displaystyle$ $\displaystyle \approx$ $\displaystyle v_{g}dP\ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{kT}{P}\;dP \ \ \ \ \ (15)$

where in the second line we used the fact that ${v_{g}\gg v_{l}}$, since the volume per gas molecule is much greater than the volume per liquid molecule. The last line assumes that the vapour is an ideal gas, so ${v_{g}=V_{g}/N_{g}=kT/P}$.

If the liquid and gas are in diffusive equilibrium, ${\mu_{g}=\mu_{l}}$ and the vapour is saturated so that ${P=P_{v}}$, so we can integrate to get

 $\displaystyle \mu_{g}-\mu_{l}$ $\displaystyle =$ $\displaystyle \int_{0}^{\mu_{g}-\mu_{l}}d\left(\mu_{g}-\mu_{l}\right)^{\prime}\ \ \ \ \ (16)$ $\displaystyle$ $\displaystyle =$ $\displaystyle kT\int_{P_{v}}^{P}\frac{dP^{\prime}}{P^{\prime}}\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle kT\ln\frac{P}{P_{v}}\ \ \ \ \ (18)$ $\displaystyle$ $\displaystyle =$ $\displaystyle kT\ln h \ \ \ \ \ (19)$

where ${h}$ is the relative humidity. Plugging this into 7 we get

$\displaystyle r_{c}=\frac{2\sigma v_{l}}{kT\ln h} \ \ \ \ \ (20)$

The snag with this equation is that the relative humidity is usually between 0 and 1, so ${\ln h<0}$ and ${r_{c}<0}$, which is clearly nonsense. In order for this equation to give a physical radius, we need ${h>1}$, so the air is ‘supersaturated’. In this case, a plot of ${r_{c}}$ versus ${h}$ looks like this:

In order for droplets to form, the initial droplet size has to be very large (relatively speaking; a radius of a few nanometres still contains a very large number of water molecules) so it would seem unlikely that such droplets would form spontaneously.

# Vapour pressure

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problem 5.41.

Two systems in diffusive equilibrium have equal chemical potentials. We can use this fact to solve the following problem.

We begin with a closed system consisting of a liquid such as water in diffusive equilibrium with its vapour. At the start, only the liquid and its vapour are present. Then we pump in an inert gas (that is, a gas that doesn’t react chemically with the liquid or vapour) to increase the pressure in the container (we’re also assuming that everything is at the same temperature, so no heat flow occurs). Assuming that the inert gas doesn’t dissolve in the liquid, and that the liquid and its vapour remain in diffusive equilibrium, the chemical potentials of the liquid and vapour must change by the same amount: ${d\mu_{\ell}=d\mu_{v}}$. What happens to the partial pressure ${P_{v}}$ of the vapour?

For the vapour, which we’ll assume is an ideal gas, the chemical potential is given by

$\displaystyle \mu_{v}\left(T,P_{v}\right)=\mu_{v}^{\circ}\left(T\right)+kT\ln\frac{P_{v}}{P^{\circ}} \ \ \ \ \ (1)$

where the superscript ‘o’ indicates the value at a reference pressure, usually taken to be 1 bar. If the pressure changes by ${dP_{v}}$, we have

$\displaystyle d\mu_{v}=\frac{kT}{P_{v}}dP_{v} \ \ \ \ \ (2)$

For the liquid, we can use the relation between the Gibbs energy and chemical potential

 $\displaystyle G$ $\displaystyle =$ $\displaystyle N\mu_{\ell} \ \ \ \ \ (3)$

If we combine this with the derivative

$\displaystyle \left(\frac{\partial G}{\partial P}\right)_{T,N}=V \ \ \ \ \ (4)$

we get

 $\displaystyle dG$ $\displaystyle =$ $\displaystyle Nd\mu_{\ell}=V\;dP\ \ \ \ \ (5)$ $\displaystyle d\mu_{\ell}$ $\displaystyle =$ $\displaystyle \frac{V}{N}dP \ \ \ \ \ (6)$

where the pressure ${P}$ is now the pressure on the liquid, which is the sum of the vapour pressure ${P_{v}}$ and the pressure due to the added inert gas. (Note that ${V}$ and ${N}$ refer to the liquid, not the vapour.) Equating 2 and 6, we get a differential equation relating the vapour pressure and total pressure:

 $\displaystyle \frac{kT}{P_{v}}dP_{v}$ $\displaystyle =$ $\displaystyle \frac{V}{N}dP\ \ \ \ \ (7)$ $\displaystyle \frac{dP_{v}}{dP}$ $\displaystyle =$ $\displaystyle \frac{V}{NkT}P_{v} \ \ \ \ \ (8)$

We can solve this as follows:

 $\displaystyle \frac{dP_{v}}{P_{v}}$ $\displaystyle =$ $\displaystyle \frac{V}{NkT}dP\ \ \ \ \ (9)$ $\displaystyle \ln P_{v}+A$ $\displaystyle =$ $\displaystyle \frac{V}{NkT}P+B \ \ \ \ \ (10)$

where ${A}$ and ${B}$ are constants of integration. To determine these constants, we use the fact that when there is no inert gas present, ${P=P_{v}}$. We’ll call this pressure ${P_{v0}}$ to distinguish it from ${P_{v}}$, the latter of which is actually a function of ${P}$. [Schroeder uses the confusing notation ${P_{v}\left(P_{v}\right)}$ to indicate the vapour pressure when there is no inert gas present, which makes it look like a function of a function.] That is, we must have

$\displaystyle \ln P_{v0}+A=\frac{V}{NkT}P_{v0}+B \ \ \ \ \ (11)$

We can choose ${A}$ and ${B}$ so that both sides are zero, giving

 $\displaystyle A$ $\displaystyle =$ $\displaystyle -\ln P_{v0}\ \ \ \ \ (12)$ $\displaystyle B$ $\displaystyle =$ $\displaystyle -\frac{V}{NkT}P_{v0} \ \ \ \ \ (13)$

Substituting these back into 10 we get

 $\displaystyle \ln\frac{P_{v}}{P_{v0}}$ $\displaystyle =$ $\displaystyle \frac{V}{NkT}\left(P-P_{v0}\right)\ \ \ \ \ (14)$ $\displaystyle P_{v}\left(P\right)$ $\displaystyle =$ $\displaystyle P_{v0}e^{\left(P-P_{v0}\right)V/NkT} \ \ \ \ \ (15)$

Thus if ${P>P_{v0}}$, the vapour pressure increases.

As an example, the vapour pressure of water at ${T=25^{\circ}\mbox{ C}}$ is ${3.169\times10^{3}\mbox{ Pa}}$. Suppose we introduce an inert gas at atmospheric pressure (${=1.01\times10^{5}\mbox{ Pa}}$) into the container. For one mole of water, ${V=18.068\times10^{-6}\mbox{ m}^{3}}$, so we have for the new vapour pressure

$\displaystyle P_{v}\left(1.01\times10^{5}\right)=3.169\times10^{3}e^{\left(1.01\times10^{5}-3.169\times10^{3}\right)\left(18.068\times10^{-6}\right)/NkT} \ \ \ \ \ (16)$

Using ${N=6.02\times10^{23}}$, ${k=1.38\times10^{-23}}$ and ${T=298}$, we have

$\displaystyle P_{v}\left(1.01\times10^{5}\right)=3.169\times10^{3}e^{7.14\times10^{-4}} \ \ \ \ \ (17)$

The exponential is very close to 1, so we can expand in a Taylor series:

$\displaystyle e^{7.14\times10^{-4}}\approx1+7.14\times10^{-4} \ \ \ \ \ (18)$

Thus the change in vapour pressure is only about 0.07%.

# Albite, jadeite and quartz phase diagram

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problem 5.40.

Here’s an example of applying the Gibbs energy and the Clausius-Clapeyron relation to a case where one solid (albite, with formula ${\mbox{NaAlSi}_{3}\mbox{O}_{8}}$) dissociates into two other solids (jadeite, with formula ${\mbox{NaAlSi}_{2}\mbox{O}_{6}}$, and quartz, ${\mbox{SiO}_{2}}$) under the reversible reaction

$\displaystyle \mbox{NaAlSi}_{3}\mbox{O}_{8}\longleftrightarrow\mbox{NaAlSi}_{2}\mbox{O}_{6}+\mbox{SiO}_{2} \ \ \ \ \ (1)$

To get a rough phase diagram, we need to first find one point on the phase boundary by using the Gibbs free energy, as we did for aluminum silicate, and then find the slope of the phase boundary using the Clausius-Clapeyron relation.

First, we use

$\displaystyle \left(\frac{\partial G}{\partial P}\right)_{T,N}=V \ \ \ \ \ (2)$

Ignoring compressibility, we can assume that ${V}$ is constant as the pressure is increased. The molar volumes, in ${\mbox{m}^{3}}$, for the three forms are, from the data in Schroeder’s book: ${V_{a}=100.07\times10^{-6}}$, ${V_{j}=60.4\times10^{-6}}$ and ${V_{q}=22.69\times10^{-6}}$. Since the pressures involved are in the kilobar region, we can take the Gibbs data (given for ${P=1\mbox{ bar}}$) to be the values for ${P\approx0}$ and we get, with energy in kJ and pressure in kbar (with ${1\mbox{ kbar}=10^{8}\mbox{ Pa}}$):

 $\displaystyle G_{a}$ $\displaystyle =$ $\displaystyle -3711.5+10.007P\ \ \ \ \ (3)$ $\displaystyle G_{j}$ $\displaystyle =$ $\displaystyle -2852.1+6.04P\ \ \ \ \ (4)$ $\displaystyle G_{q}$ $\displaystyle =$ $\displaystyle -856.4+2.269P \ \ \ \ \ (5)$

For the reaction 1, the relevant Gibbs energy for the RHS is the sum ${G_{j}+G_{q}}$, so the energies we want to compare to find the stable form are

 $\displaystyle G_{a}$ $\displaystyle =$ $\displaystyle -3711.5+10.007P\ \ \ \ \ (6)$ $\displaystyle G_{jq}$ $\displaystyle =$ $\displaystyle -3708.5+8.309P \ \ \ \ \ (7)$

Plotting these two lines we get

The red line is albite and the green line is jadeite + quartz. At low pressure, albite has the lower Gibbs energy so is the stable form. The lines cross at a pressure of ${P_{ajq}=1.77\mbox{ kbar}}$ so at that pressure and a temperature of 298 K, the two sides are in equilibrium. For higher pressures, the stable form is jadeite + quartz.

We could use this point as one point on the phase diagram, but we can alse use the earlier method that we applied to aluminum silicate. The difference in Gibbs energy between the two sides of the equation is

$\displaystyle \Delta G\left(T_{2}\right)=\Delta G\left(T_{1}\right)-\Delta S\left(T_{2}-T_{1}\right) \ \ \ \ \ (8)$

The entropies are, in units of ${\mbox{kJ K}^{-1}}$: ${S_{a}=0.2074}$, ${S_{j}=0.1335}$ and ${S_{q}=0.04184}$ so, taking ${T_{1}=298\mbox{ K}}$ and ${\Delta G\left(298\right)=G_{a}-G_{jq}}$ evaluated at ${P=0}$ from 6, we have

$\displaystyle \Delta G\left(T_{2}\right)=-3.0-0.03206\left(T-298\right) \ \ \ \ \ (9)$

The phase boundary passes through the point where ${\Delta G=0}$, which occurs at ${T=204.4\mbox{ K}}$. From the Clausius-Clapeyron equation, the slope of the boundary at this point is

$\displaystyle \frac{dP}{dT}=\frac{S_{a}-\left(S_{j}+S_{q}\right)}{V_{a}-\left(V_{j}+V_{q}\right)}=0.0189\mbox{ kbar K}^{-1} \ \ \ \ \ (10)$

Thus the phase boundary has the equation

 $\displaystyle P$ $\displaystyle =$ $\displaystyle \frac{dP}{dT}\left(T-204.4\right)\ \ \ \ \ (11)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 0.0189\left(T-204.4\right) \ \ \ \ \ (12)$

A plot of the boundary looks like this:

As always, this is only a crude approximation as it disregards variations in entropy and volume as the pressure and temperature are changed, and extrapolates the entire boundary from the slope at a single point.

# Phases of water – plots of Gibbs energy

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problems 5.30 – 5.31.

We can apply a similar analysis to that for graphite and diamond to the three phases of water.

Using the data in Schroeder’s book, we can get some insight into the stability of these three forms by considering the Gibbs free energy. The thermodynamic potential for the Gibbs free energy is

$\displaystyle dG=-S\;dT+V\;dP+\mu\;dN \ \ \ \ \ (1)$

If ${P}$ and ${N}$ are constant, we have

$\displaystyle \left(\frac{\partial G}{\partial T}\right)_{P,N}=-S \ \ \ \ \ (2)$

If ${T}$ and ${N}$ are constant, we have

$\displaystyle \left(\frac{\partial G}{\partial P}\right)_{T,N}=V \ \ \ \ \ (3)$

We can use these derivatives to sketch qualitative graphs for ${G}$ as a function of ${P}$ and ${T}$ for water. First, we look at 2. Since the slope of the ${G}$ versus ${T}$ curve is the negative entropy, and the entropy of steam is higher than liquid water which is in turn higher than ice, we’d expect the curves for ice, water and steam to have increasingly negative slopes. To arrange them in the correct relative position, we can look at a phase diagram for water, such as Schroeder’s Figure 5.11. At a pressure of ${P=1\mbox{ bar}}$, ice and liquid water are equally stable at ${T=0^{\circ}\mbox{ C}}$, so we’d expect the ${G}$ curves for water and ice to intersect at that temperature. Below ${0^{\circ}\mbox{ C}}$, ice is the stable form. At ${T=100^{\circ}\mbox{ C}}$, water and steam are equally stable, so we’d expect the ${G}$ curves for water and ice to intersect at that temperature, with steam being the stable form above ${T=100^{\circ}\mbox{ C}}$. Thus the plot should look something like this (in all these plots, the values on the ${G}$ axis are not numerically correct; they merely serve to show whether ${G}$ is higher or lower when comparing curves):

Here, red is ice, green is liquid water and blue is steam.

For a much smaller pressure such as ${P=0.001\mbox{ bar}}$, the phase diagram shows that ice sublimes directly to steam at a temperature less than ${0^{\circ}\mbox{ C}}$ so the graph would look something like this:

At this pressure, water is never found in liquid form.

Turning now to 3, we see that at constant temperature, the slope of the ${G}$ versus ${P}$ plot gives the volume. Consider first ${T=0^{\circ}\mbox{ C}}$. At low pressure, we have only steam which has a large molar volume, so we’d expect the slope to be large and positive (tending to ${+\infty}$ as ${P\rightarrow0}$). From the phase diagram, we see that as we increase the pressure, the steam first condenses to ice, then at ${P=1\mbox{ bar}}$ (roughly), ice and water are equally stable. Increasing the pressure further actually causes the ice to melt into liquid water (a peculiarity of water, since its liquid phase is denser than the solid; with most substances, increasing the pressure converts a liquid into a solid). Thus we’d expect a plot something like this:

There’s a small interval where the solid phase (red) is the most stable.

If we increase the temperature a few degrees, then we still have steam at low pressures but this condenses directly to liquid as we increase the pressure, so there is no point where the solid is the stable form. The graph looks like this:

We see that the solid (red line) is never the lowest energy, so is never stable.

# Aluminum silicate – stability of three crystal structures

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problem 5.29.

We can apply a similar analysis to that for graphite and diamond to the more complex case of aluminum silicate ${\mbox{Al}_{2}\mbox{SiO}_{3}}$, which has three crystalline forms: kyanite, andalusite and sillimanite. Each is stable under some range of temperature and pressure.

Using the data in Schroeder’s book, we can get some insight into the stability of these three forms by considering the Gibbs free energy. The thermodynamic potential for the Gibbs free energy is

$\displaystyle dG=-S\;dT+V\;dP+\mu\;dN \ \ \ \ \ (1)$

First, if ${T}$ and ${N}$ are constant, we have

$\displaystyle \left(\frac{\partial G}{\partial P}\right)_{T,N}=V \ \ \ \ \ (2)$

Ignoring compressibility, we can assume that ${V}$ is constant as the pressure is increased. The molar volumes for the three forms are ${V_{k}=44.09\times10^{-6}\mbox{ m}^{3}}$ for kyanite, ${V_{a}=51.53\times10^{-6}\mbox{ m}^{3}}$ for andalusite and ${V_{s}=49.90\times10^{-6}\mbox{ m}^{3}}$ for sillimanite. Therefore the Gibbs energies are

 $\displaystyle G_{k}$ $\displaystyle =$ $\displaystyle G_{k0}+44.09\times10^{-6}P\ \ \ \ \ (3)$ $\displaystyle G_{a}$ $\displaystyle =$ $\displaystyle G_{a0}+51.53\times10^{-6}P\ \ \ \ \ (4)$ $\displaystyle G_{s}$ $\displaystyle =$ $\displaystyle G_{s0}+49.90\times10^{-6}P \ \ \ \ \ (5)$

The constants of integration ${G_{k0}}$, ${G_{a0}}$ and ${G_{s0}}$ can be found by requiring the energies to be those quoted in Schroeder for ${P=1\mbox{ bar}=10^{5}\mbox{ Pa}}$. As with our calculations for calcium carbonate, the difference in Gibbs energy between a pressure of zero and 1 bar is negligible so we can just set the constants of integration to the values quoted in Schroeder for ${P=1\mbox{ bar}}$. We get

 $\displaystyle G_{k}$ $\displaystyle =$ $\displaystyle -2443.88\times10^{3}+44.09\times10^{-6}P\ \ \ \ \ (6)$ $\displaystyle G_{a}$ $\displaystyle =$ $\displaystyle -2442.66\times10^{3}+51.53\times10^{-6}P\ \ \ \ \ (7)$ $\displaystyle G_{s}$ $\displaystyle =$ $\displaystyle -2440.99\times10^{3}+49.90\times10^{-6}P \ \ \ \ \ (8)$

At ${P=0}$, ${G_{k}}$ is the lowest and since the slope of the line in the plot of ${G_{k}}$ versus ${P}$ is also the smallest, ${G_{k}}$ will be less than both ${G_{a}}$ and ${G_{s}}$ at all pressures, so at a temperature of 298 K, kyanite is always the most stable form. The graph looks like this:

The red line is kyanite, green is andalusite and yellow is sillimanite. The latter two lines do intersect at a pressure of around 9 kbar, but kyanite is always the most stable.

At constant pressure, we can use 1 to investigate the stability of the three forms as temperature varies. We have

$\displaystyle \left(\frac{\partial G}{\partial T}\right)_{P,N}=-S \ \ \ \ \ (9)$

By subtracting the two forms of this equation for two crystal types, we can get a relation for the difference in Gibbs energy between two types:

$\displaystyle \left(\frac{\partial\Delta G}{\partial T}\right)_{P,N}=-\Delta S \ \ \ \ \ (10)$

The entropy of a substance can be written as

$\displaystyle dS=\frac{dQ}{T}=C_{P}\frac{dT}{T} \ \ \ \ \ (11)$

so the total entropy as a function of temperature depends on the heat capacity ${C_{P}}$. If we take Schroeder’s word for it that ${C_{P}}$ at high temperatures depends only on ${N}$, then even though the raw value of ${S}$ for a sample of solid increases with ${T}$, the difference between the entropy values for two substances should remain relatively constant as we increase ${T}$. We can therefore integrate 10 to get

 $\displaystyle \Delta G\left(T_{2}\right)$ $\displaystyle =$ $\displaystyle \Delta G\left(T_{1}\right)-\int_{T_{1}}^{T_{2}}\Delta S\left(T\right)\;dT\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \Delta G\left(T_{1}\right)-\Delta S\left(T_{2}-T_{1}\right) \ \ \ \ \ (13)$

The entropy values quoted in Schroeder’s book for 1 mole at ${T=298\mbox{ K}}$ and ${P=1\mbox{ bar}}$ are ${S_{k}=83.81\mbox{ J K}^{-1}}$, ${S_{a}=93.22\mbox{ J K}^{-1}}$ and ${S_{s}=96.11\mbox{ J K}^{-1}}$ so using the Gibbs energy values from above, we get (with energies in kJ):

 $\displaystyle \Delta G_{ka}$ $\displaystyle =$ $\displaystyle G_{k}-G_{a}-\left(S_{k}-S_{a}\right)\left(T-298\right)\ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -1.22+9.41\times10^{-3}\left(T-298\right)\ \ \ \ \ (15)$ $\displaystyle \Delta G_{ks}$ $\displaystyle =$ $\displaystyle -2.89+12.30\times10^{-3}\left(T-298\right)\ \ \ \ \ (16)$ $\displaystyle \Delta G_{as}$ $\displaystyle =$ $\displaystyle \Delta G_{ks}-\Delta G_{ka}\ \ \ \ \ (17)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -1.67+2.89\times10^{-3}\left(T-298\right) \ \ \ \ \ (18)$

If we plot these three lines, we get

Here, ${\Delta G_{ka}}$ is red, ${\Delta G_{ks}}$ is green and ${\Delta G_{as}}$ is yellow. To use this graph to determine the temperature range for which each structure is the most stable, we observe that the most stable form has the lowest value of ${G}$ at that temperature. For kyanite to be the most stable, for example, we must have ${G_{k} and ${G_{k}. Therefore, ${\Delta G_{ka}=G_{k}-G_{a}<0}$ and ${\Delta G_{ks}=G_{k}-G_{s}<0}$. Looking at the graph, these conditions are true if both the red and green lines are below the ${T}$ axis, which occurs when ${T}$ is less than the point where the red line crosses the ${T}$ axis, which is where ${\Delta G_{ka}=0}$. From 15 this occurs at ${T=427.6\mbox{ K}}$.

For andalusite to be stable, we must have ${G_{a} and ${G_{a}, which means ${\Delta G_{ka}>0}$ and ${\Delta G_{as}<0}$. Thus the red line must be above the ${T}$ axis and the yellow line must be below the ${T}$ axis. The minimum ${T}$ satisfying these conditions is ${T=427.6\mbox{ K}}$ that we found above. The maximum ${T}$ is where the yellow line crosses the ${T}$ axis, which occurs when ${\Delta G_{as}=0}$, at ${T=875.9\mbox{ K}}$.

For sillimanite to be stable, we must have ${G_{s} and ${G_{s}, which means ${\Delta G_{ks}>0}$ and ${\Delta G_{as}>0}$. Thus both the green and yellow lines must be above the ${T}$ axis, which is true for ${T>875.9\mbox{ K}}$. To summarize: kyanite is stable for ${T<427.6\mbox{ K}}$, andalusite is stable for ${427.6\mbox{ K} and sillimanite is stable for ${T>875.9\mbox{ K}}$.

How accurate is the assumption that ${\Delta S}$ (the difference in entropy between two states of the substance) is actually a constant over a wide range of temperature? From 11 for one state, we can integrate:

 $\displaystyle S\left(T_{2}\right)$ $\displaystyle =$ $\displaystyle S\left(T_{1}\right)+C_{P}\int_{T_{1}}^{T_{2}}\frac{dT}{T}\ \ \ \ \ (19)$ $\displaystyle$ $\displaystyle =$ $\displaystyle S\left(T_{1}\right)+C_{P}\ln\frac{T_{2}}{T_{1}} \ \ \ \ \ (20)$

Thus ${\Delta S}$ between two states, ${i}$ and ${j}$ is obtained by subtracting one equation of this type from the other:

$\displaystyle \Delta S_{ij}\left(T_{2}\right)=\Delta S_{ij}\left(T_{1}\right)+\left(C_{Pi}-C_{Pj}\right)\ln\frac{T_{2}}{T_{1}} \ \ \ \ \ (21)$

${\Delta S_{ij}}$ will be constant if the two states have equal heat capacities. From the data in Schroeder’s book, ${C_{Pk}=121.71\mbox{ J K}^{-1}}$, ${C_{Pa}=122.72\mbox{ J K}^{-1}}$ and ${C_{Ps}=124.52\mbox{ J K}^{-1}}$. Using a reference temperature of ${T_{1}=298\mbox{ K}}$, the log factor over the range of temperatures over which the state changes occur varies from ${\ln\frac{427.6}{298}=0.361}$ to ${\ln\frac{875.9}{298}=1.08}$. The largest difference in heat capacity is ${C_{Ps}-C_{Pk}=2.81}$, so we can expect the last term in 21 to be in the vicinity of ${1\mbox{ J K}^{-1}}$. The entropy differences used in 15 and following equations are between around 2 and 12 ${\mbox{J K}^{-1}}$ so actual value of ${\Delta S}$ could be in error by a fair bit. Of course, we’ve also used the heat capacities for ${T=298\mbox{ K}}$ and assumed that they don’t vary over the temperature range, so we’d really need to correct that error as well.

# Calcite and aragonite

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problem 5.28.

We can apply a similar analysis to that for graphite and diamond to calcium carbonate ${\mbox{CaCO}_{3}}$, which also has two crystalline forms, calcite and aragonite. Calcite is a clear, colourless crystal which has the curious property of double refraction. Here’s an image of a calcite crystal placed on Schroeder’s problem 5.28:

You can see that a double image is produced.

Using the data in Schroeder’s book, we can get some insight into the stability of these two forms by considering the Gibbs free energy. From the thermodynamic potential for the Gibbs free energy

$\displaystyle dG=-S\;dT+V\;dP+\mu\;dN \ \ \ \ \ (1)$

we have

$\displaystyle \left(\frac{\partial G}{\partial P}\right)_{T,N}=V \ \ \ \ \ (2)$

Ignoring compressibility, we can assume that ${V}$ is constant as the pressure is increased. The molar volumes are ${V_{c}=36.93\times10^{-6}\mbox{ m}^{3}}$ for calcite and ${V_{a}=34.15\times10^{-6}\mbox{ m}^{3}}$ for aragonite. Therefore the Gibbs energies are

 $\displaystyle G_{c}$ $\displaystyle =$ $\displaystyle G_{c0}+36.93\times10^{-6}P\ \ \ \ \ (3)$ $\displaystyle G_{a}$ $\displaystyle =$ $\displaystyle G_{a0}+34.15\times10^{-6}P \ \ \ \ \ (4)$

The constants of integration ${G_{c0}}$ and ${G_{a0}}$ can be found by requiring the energies to be those quoted in Schroeder for ${P=1\mbox{ bar}=10^{5}\mbox{ Pa}}$. The energies are ${-1128.8\times10^{3}\mbox{ J}}$ for calcite and ${-1127.8\times10^{3}\mbox{ J}}$ for aragonite. With these values, we have

 $\displaystyle G_{c0}$ $\displaystyle =$ $\displaystyle -1128.8\times10^{3}-3.693\approx-1128.8\times10^{3}\mbox{ J}\ \ \ \ \ (5)$ $\displaystyle G_{a0}$ $\displaystyle =$ $\displaystyle -1127.8\times10^{3}-3.415\approx-1127.8\times10^{3}\mbox{ J} \ \ \ \ \ (6)$

So the energies are

 $\displaystyle G_{c}$ $\displaystyle =$ $\displaystyle -1128.8\times10^{3}+36.93\times10^{-6}P\ \ \ \ \ (7)$ $\displaystyle G_{a}$ $\displaystyle =$ $\displaystyle -1127.8\times10^{3}+34.15\times10^{-6}P \ \ \ \ \ (8)$

Since ${G_{c}\left(P=10^{5}\right), calcite is (very slightly) more stable than aragonite at atmospheric pressure and room temperature. To find the pressure at which aragonite becomes more stable, we set the two energies equal to each other:

 $\displaystyle -1128.8\times10^{3}+36.93\times10^{-6}P$ $\displaystyle =$ $\displaystyle -1127.8\times10^{3}+34.15\times10^{-6}P\ \ \ \ \ (9)$ $\displaystyle P$ $\displaystyle =$ $\displaystyle 3.597\times10^{8}\mbox{ Pa}\ \ \ \ \ (10)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 3.597\mbox{ kbar} \ \ \ \ \ (11)$

Because the volumes are very similar, a plot of the two Gibbs energies shows two lines that are nearly parallel. In the plot, calcite is in red and aragonite is in blue.

# Graphite and diamond

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problems 5.24 – 5.27.

Two of the main forms of elemental carbon are graphite and diamond. We can get some insight into the stability of these two forms by considering the Gibbs free energy. From the thermodynamic potential for the Gibbs free energy

$\displaystyle dG=-S\;dT+V\;dP+\mu\;dN \ \ \ \ \ (1)$

we have

$\displaystyle \left(\frac{\partial G}{\partial P}\right)_{T,N}=V \ \ \ \ \ (2)$

Using data given in Schroeder’s book, the molar volume of graphite is ${V_{g}=5.31\times10^{-6}\mbox{ m}^{3}}$ and of diamond is ${V_{d}=3.42\times10^{-6}\mbox{ m}^{3}}$. If we ignore compression effects over a wide range of pressures, we can integrate 2 to get a linear relationshipe between ${G}$ and ${P}$. If we take the Gibbs energy of graphite to be ${G_{g}=0}$ at ${T=298\mbox{ K}}$ and ${P=10^{5}\mbox{ Pa}}$ then for diamond under the same conditions we have ${G_{d}=2900\mbox{ J}}$, so the relations are (for one mole of carbon atoms):

 $\displaystyle G_{g}$ $\displaystyle =$ $\displaystyle V_{g}P=5.31\times10^{-6}P\ \ \ \ \ (3)$ $\displaystyle G_{d}$ $\displaystyle =$ $\displaystyle V_{d}P=2900+3.42\times10^{-6}P \ \ \ \ \ (4)$

As the diamond line starts above the graphite line but has a smaller slope, it will intersect the graphite line at a pressure ${P_{gd}}$ determined by

 $\displaystyle 5.31\times10^{-6}P_{gd}$ $\displaystyle =$ $\displaystyle 2900+3.42\times10^{-6}P_{gd}\ \ \ \ \ (5)$ $\displaystyle P_{gd}$ $\displaystyle =$ $\displaystyle 1.534\times10^{9}\mbox{ Pa}\ \ \ \ \ (6)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 15.34\mbox{ kbar} \ \ \ \ \ (7)$

using the conversion ${10^{5}\mbox{ Pa}=1\mbox{ bar}}$.

If we include the effect of compressibility, things get a bit more complicated. The isothermal compressibility of graphite (given by Schroeder) is about

$\displaystyle \kappa_{T}=-\frac{1}{V}\left(\frac{\partial V}{\partial P}\right)_{T,N}=3\times10^{-6}\mbox{ bar}^{-1}=3\times10^{-11}\mbox{ Pa}^{-1} \ \ \ \ \ (8)$

Assuming ${\kappa_{T}}$ is constant over the range of pressures we’re considering, we can integrate this to get

 $\displaystyle \frac{dV}{V}$ $\displaystyle =$ $\displaystyle -\kappa_{T}dP\ \ \ \ \ (9)$ $\displaystyle \ln\frac{V}{V_{0}}$ $\displaystyle =$ $\displaystyle \kappa_{T}\left(P_{0}-P\right) \ \ \ \ \ (10)$

where ${V_{0}}$ and ${P_{0}}$ are some reference volume and pressure, which we can take as the molar volume and a pressure of 1 bar.

Thus the volume is a function of the pressure:

$\displaystyle V=V_{0}e^{\kappa_{T}P_{0}}e^{-\kappa_{T}P} \ \ \ \ \ (11)$

For graphite, this gives

 $\displaystyle V_{g}$ $\displaystyle =$ $\displaystyle \left(5.31\times10^{-6}\right)e^{3\times10^{-6}}e^{-\kappa_{T}P}\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle \approx$ $\displaystyle 5.31\times10^{-6}e^{-3\times10^{-11}P} \ \ \ \ \ (13)$

since ${e^{3\times10^{-6}}\approx1}$. The Gibbs energy for graphite can then be found from 2 as

 $\displaystyle G_{g}$ $\displaystyle =$ $\displaystyle -\frac{5.31\times10^{-6}}{3\times10^{-11}}e^{-3\times10^{-11}P}+G_{0}\ \ \ \ \ (14)$ $\displaystyle$ $\displaystyle =$ $\displaystyle -1.77\times10^{5}e^{-3\times10^{-11}P}+G_{0} \ \ \ \ \ (15)$

The integration constant ${G_{0}}$ can be found from the condition that ${G_{g}\left(10^{5}\mbox{ Pa}\right)=0}$, so ${G_{0}=1.77\times10^{5}\mbox{ J}}$ and we get

$\displaystyle G_{g}=1.77\times10^{5}\left(1-e^{-3\times10^{-11}P}\right) \ \ \ \ \ (16)$

As ${\kappa_{T}}$ for diamond is about ten times smaller than ${\kappa_{T}}$ for graphite, we can approximate ${G_{d}}$ by the same straight line that we used above. The pressure at which diamond becomes more stable than graphite is found by equating ${G_{g}}$ and ${G_{d}}$ as before, although this time we get a transcendental equation (with ${P}$ found both in an exponent and as a linear factor), so we can solve the equation numerically using Maple:

 $\displaystyle 1.77\times10^{5}\left(1-e^{-3\times10^{-11}P_{gd}}\right)$ $\displaystyle =$ $\displaystyle 2900+3.42\times10^{-6}P_{gd}\ \ \ \ \ (17)$ $\displaystyle P_{gd}$ $\displaystyle =$ $\displaystyle 1.647\times10^{9}\mbox{ Pa}\ \ \ \ \ (18)$ $\displaystyle$ $\displaystyle =$ $\displaystyle 16.47\mbox{ kbar} \ \ \ \ \ (19)$

The Gibbs energies look like this, with the blue curve being graphite and the red curve being diamond:

The blue curve is actually slightly concanve downwards, as can be seen by plotting it for higher pressures, but in the pressure range here, it’s approximately linear. The red curve is actually a straight line.

There are a couple of rather odd features about the graphite-diamond transition. First, from the Gibbs energy plot, we see that graphite is actually the more stable form at atmospheric pressure and temperature, so that should mean that diamonds will spontaneously degrade into graphite over time. Fortunately for owners of diamonds, this process takes a very long time (although it should also be noted that diamonds, being carbon, do burn, so don’t throw your wedding ring into a fire!).

Another odd fact is that because of its more ordered crystal structure, diamond has a lower entropy than graphite (for 1 mole, the values are ${S_{d}=2.38\mbox{ J K}^{-1}}$ and ${S_{g}=5.74\mbox{ J K}^{-1}}$). Schroeder asks us to explain how the conversion of graphite to diamond at high pressure can result in an overall increase of the environment’s entropy. It’s not clear to me why this should be, though perhaps because ${G_{d} at high pressure, the energy transferred away from the graphite into the environment causes an overall increase in entropy. Comments welcome.

Finally, geochemists use a bizarre unit of volume equal to ${\mbox{kJ kbar}^{-1}}$. The conversion to the more standard cubic metre is given by

$\displaystyle 1\mbox{ kJ kbar}^{-1}=1\mbox{ J bar}^{-1}=10^{-5}\mbox{ J Pa}^{-1}=10^{-5}\mbox{ m}^{3} \ \ \ \ \ (20)$