Linear compressibility Lennard-Jones potential

[WannabeNewton]

Homework Statement

Please see attachment. The correct answer is supposed to be ##\alpha = .15 k_B/U_0##.

The Attempt at a Solution

The minimum of the potential is at ##x = a##. Expanded about this minimum we get ##U = -U_0 + 36\frac{U_0}{a^2}(x - a)^2 - 252\frac{U_0}{a^3}(x - a)^3 + ...##. The mean distance is ##\langle x \rangle = \frac{\int_{0}^{\infty} x e^{-\beta U} dx}{\int_{0}^{\infty}e^{-\beta U}dx}##. Here is a plot of the potential: http://en.wikipedia.org/wiki/Lenard_Jones_potential and as you can see it falls of steeply from infinity at ##x =0## to the minimum at ##x = a## and then rises slowly to zero as ##x\rightarrow \infty##. In order to insert the above series expansion of ##U## about ##x = a## into the expression for ##\langle x \rangle##, we would need to be able to show that the integrals in it are overwhelmingly dominated by the values of the integrand near ##x = a##; let's say that there is a characteristic width ##\delta a## that codifies being near ##x = a##.

As we move away from ##x = a## and towards ##x = 0##, we see that ##U/U_0## very, very rapidly becomes larger than unity so we can say to good approximation that ##U/U_0 \gg 1## when ##x \leq a - \delta a##. Furthermore ##\beta U_0 \gg 1## by hypothesis so ##\beta U \gg 1## and note ##U > 0## when ##x \leq a - \delta a##. Hence ##\int_{0}^{a - \delta a} x e^{-\beta U} dx \approx 0, \int_{0}^{a - \delta a} e^{-\beta U} dx \approx 0##.

Now we have to show that ##\int_{a + \delta a}^{\infty} x e^{-\beta U}dx \approx 0, \int_{a+ \delta a}^{\infty}e^{-\beta U}dx \approx 0##. This is where I'm not sure. We see that for ##x \geq a + \delta a##, ##U/U_0 \ll 1##, ##U < 0##, and ##U \rightarrow 0##. So we see that even if ##\beta U_0 \gg 1##, eventually ##\beta U \ll 1## i.e. for ##x \gg a + \delta a##, which from the diagram happens as soon as ##x \approx 2##, but ##U < 0## so ##e^{-\beta U} = e^{|\beta U|} \approx 0## for ##x \geq a + \delta a## and what's more ##e^{|\beta U|} \rightarrow 0## much more rapidly than ##x \rightarrow \infty## so we get the desired results.
Thus ##\langle x \rangle = \frac{\int_{0}^{\infty} x e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} e^{252\beta U_0 \frac{(x - a)^3}{a^3}}...dx}{\int_{0}^{\infty} e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} e^{252\beta U_0 \frac{(x - a)^3}{a^3}}...dx}##.

We know that ##\delta a/a \ll 1## so it would be natural to try working only to second order in ##\frac{(x - a)}{a}## since ##x - a \sim \delta a## in the integrals. To second order the integrals are trivial to evaluate but give the wrong result of ##\alpha = .5/T##. I don't know why. The next step then, albeit a very ad-hoc step, would be to try working to third order but beyond second order, the order at which we work for this particular problem seems quite arbitrary rather than natural. At this point I'm just trying the third order in order to get what the book says is the right answer. I don't have any physical justification for why the third order would be the right order to work in, as opposed to second order (which for some reason gives the wrong answer) or any higher orders.

Put in other words, for this problem how would I go about figuring out to what order in ##\frac{(x - a)}{a}## I should work in? I don't want to just keep trying higher orders until I get the right answer, there has to be a more physically intuitive way to figure it out. Could anyone shed some light on this issue? Does it have something to do with the size of ##\pm \beta U_0 \frac{(x - a)}{a} \sim \pm \beta U_0 \frac{\delta a}{a}## as opposed to just the size of ##\frac{(x - a)}{a}## itself? Thanks in advance.

P.S. I avoided posting the third order calculation in the OP itself because then the OP would just get way too long to read.

 

[vanhees71]

To answer the question about the validity of the approximation, you can use the "saddle-point approximation" for the integral over the exponential function. It's a pretty useful technique. Unfortunately, I don't know a good English book about it, and my own manuscript on "math for physicists" is in German either.

It's treated somewhat in

A. Sommerfeld, Lectures on Theoretical Physics, Vol. 6 (Partial differential equations)

That's anyway a marvelous book on mathematical methods in theoretical physics, covering Fourier and generalized Fourier series and a lot of nice tricks using complex-function theory.

 

[TSny]

Approximating only to the quadratic terms will not give any thermal expansion. A quadratic potential well is symmetric about the minimum. No matter what the amplitude of the oscillations about the minimum, the average value of x will not change. So, you need to include the cubic term.

The saddle point method also occurred to me. I'm not proficient at using it. I did find a treatment here. I tried using the approximation as given in (1.7) along with the expression for C₂ given below (1.7). Also see (1.9). I could not get the right numerical coefficient of .15 in the answer.

 

[WannabeNewton]

Thanks for the replies guys.

No matter what the amplitude of the oscillations about the minimum, the average value of x will not change.

This would imply that ##\frac{\partial}{\partial T}\langle x \rangle = 0## to second order which is not true in the case at hand. To second order it comes out to ##\frac{\partial}{\partial T}\langle x \rangle = .5/T##. If it vanished then yes I agree the third order is necessary but it doesn't so I'm afraid I'm not following you.

The saddle point method also occurred to me. I'm not proficient at using it. I did find a treatment here. I tried using the approximation as given in (1.7) along with the expression for C₂ given below (1.7). Also see (1.9). I could not get the right numerical coefficient of .15 in the answer.

Alright so starting from the very top we have the two integrals ##\int_{0}^{\infty} e^{\beta U_0 g(x)}dx## and ##\int_{0}^{\infty}x e^{\beta U_0 g(x)}dx## to be evaluated, where ##g(x) = -[(a/x)^{12} - 2(a/x)^{6}]##. We would therefore have to apply the saddle point method twice, once with ##f(x) = 1## and once with ##f(x) = x## using ##A = \beta U_0 \gg 1##. The maximum of ##g(x)## is just the minimum of the potential modulo the factor of ##U_0## so the maximum occurs at ##x_0 = a## with ##g(x_0) = 1##.

Applying (1.7) for ##f(x) = 1## to lowest non-vanishing order in ##\beta U_0## gives ##e^{\beta U_0}\sqrt{\frac{\pi a^2}{36 \beta U_0}}(1 + \frac{C_{2}}{\beta U_0}) + O(\beta^{-2} U_0^{-2})## where ##C_2 = g''''/(32 (g'')^2) - 5(g''')^2/(192(g'')^3)## evaluated at ##x_0 = a## so ##C_2 \approx .0015##.

For ##f(x) = x## to lowest non-vanishing order in ##\beta U_0## we have ##a e^{\beta U_0}\sqrt{\frac{\pi a^2}{36 \beta U_0}}(1 + \frac{C_{2}}{\beta U_0}) + O(\beta^{-2} U_0^{-2})## with ##C_2 = g'''/(8a_0 (g'')^2) + .0015 \approx .038##

So finally ##\langle x \rangle = a \frac{1 + \frac{.038}{\beta U_0}}{1 + \frac{.0015}{\beta U_0}}##. Expanding again to lowest non-vanishing order in ##\beta U_0## gives ##\langle x \rangle \approx a_0 (1 + \frac{.036}{\beta U_0})## yielding ##\alpha = \frac{1}{\langle x \rangle}\frac{\partial \langle x \rangle}{\partial T} = \frac{.036 k_B/U_0}{1 + \frac{.036}{\beta U_0}} = .036k_B/U_0 + O(\beta^{-2} U_0^{-2})##. Is this what you got?

 

[TSny]

This would imply that ##\frac{\partial}{\partial T}\langle x \rangle = 0## to second order which is not true in the case at hand. To second order it comes out to ##\frac{\partial}{\partial T}\langle x \rangle = .5/T##. If it vanished then yes I agree the third order is necessary but it doesn't so I'm afraid I'm not following you.

I don't see how you got this result.
If I keep only the quadratic approximation and use the fact that U_o >> kT, then I get virtually no T dependence of <x>. The approximation U_o >> kT allows you to extend the integration region for the quadratic approximation from (0, ∞) to (-∞, ∞).

Alright so starting from the very top we have the two integrals ##\int_{0}^{\infty} e^{\beta U_0 g(x)}dx## and ##\int_{0}^{\infty}x e^{\beta U_0 g(x)}dx## to be evaluated, where ##g(x) = -[(a/x)^{12} - 2(a/x)^{6}]##. We would therefore have to apply the saddle point method twice, once with ##f(x) = 1## and once with ##f(x) = x## using ##A = \beta U_0 \gg 1##. The maximum of ##g(x)## is just the minimum of the potential modulo the factor of ##U_0## so the maximum occurs at ##x_0 = a## with ##g(x_0) = 1##.

Applying (1.7) for ##f(x) = 1## to lowest non-vanishing order in ##\beta U_0## gives ##e^{\beta U_0}\sqrt{\frac{\pi a^2}{36 \beta U_0}}(1 + \frac{C_{2}}{\beta U_0}) + O(\beta^{-2} U_0^{-2})## where ##C_2 = g''''/(32 (g'')^2) - 5(g''')^2/(192(g'')^3)## evaluated at ##x_0 = a## so ##C_2 \approx .0015##.

For ##f(x) = x## to lowest non-vanishing order in ##\beta U_0## we have ##a e^{\beta U_0}\sqrt{\frac{\pi a^2}{36 \beta U_0}}(1 + \frac{C_{2}}{\beta U_0}) + O(\beta^{-2} U_0^{-2})## with ##C_2 = g'''/(8a_0 (g'')^2) + .0015 \approx .038##

So finally ##\langle x \rangle = a \frac{1 + \frac{.038}{\beta U_0}}{1 + \frac{.0015}{\beta U_0}}##. Expanding again to lowest non-vanishing order in ##\beta U_0## gives ##\langle x \rangle \approx a_0 (1 + \frac{.036}{\beta U_0})## yielding ##\alpha = \frac{1}{\langle x \rangle}\frac{\partial \langle x \rangle}{\partial T} = \frac{.036 k_B/U_0}{1 + \frac{.036}{\beta U_0}} = .036k_B/U_0 + O(\beta^{-2} U_0^{-2})##. Is this what you got?

That's what I get, too.

 

[TSny]

##\langle x \rangle = \frac{\int_{0}^{\infty} x e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} e^{252\beta U_0 \frac{(x - a)^3}{a^3}}...dx}{\int_{0}^{\infty} e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} e^{252\beta U_0 \frac{(x - a)^3}{a^3}}...dx}##

If I forget about the saddle point method and just approximate ##e^{\left(252 \beta U_0 (x-a)^3/a^3\right)}## using ##e^x \approx 1+x##, I get the 0.15 coefficient. The integration region may be extended to (-∞, ∞).

Don't know why the saddle point method gives a different result.

 

[WannabeNewton]

I don't see how you got this result.

Well to second order we have [tex]\langle x \rangle = \frac{\int_{0}^{\infty} x e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} dx}{\int_{0}^{\infty} e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} dx} \approx \frac{\int_{a}^{\infty} x e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} dx}{\int_{a}^{\infty} e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} dx}\\ = a + \frac{\int_{0}^{\infty} x e^{-36\beta U_0 \frac{x^2}{a^2}} dx}{\int_{0}^{\infty} e^{-36\beta U_0 \frac{x^2}{a^2}} dx} = a[1 + \frac{1}{6\sqrt{\pi}}(1/\beta U_0)^{1/2}][/tex] so that [tex]\alpha = \frac{1}{\langle x \rangle}\frac{\partial \langle x \rangle}{\partial T} = \frac{\frac{1}{12(\pi)^{3/2}}\frac{k_B}{U_0}(\beta U_0)^{1/2}}{1 + \frac{1}{6\sqrt{\pi}}(1/\beta U_0)^{1/2}} \approx \frac{1}{12(\pi)^{3/2}}\frac{k_B}{U_0}(\beta U_0)^{1/2} - \frac{1}{72 \pi^2}\frac{k_B}{U_0}[/tex] which I now realize is actually much worse in form than what I originally said I got but it certainly has a temperature dependence to leading order.

If I forget about the saddle point method and just approximate ##e^{\left(252 \beta U_0 (x-a)^3/a^3\right)}## using ##e^x \approx 1+x##, I get the 0.15 coefficient. The integration region may be extended to (-∞, ∞).

I understand why we can extend the integral to ##-\infty## to an excellent approximation but why is it valid to expand ##e^{252 \beta U_0 (x-a)^3/a^3 }## to first order? This is what I was originally stuck on actually. In the integral itself we have the product ##e^{252 \beta U_0 (x-a)^3/a^3 }e^{-36 \beta U_0 (x-a)^2/a^2 }##. So why expand the ## e^{252 \beta U_0 (x-a)^3/a^3 }## term and not the ##e^{-36 \beta U_0 (x-a)^2/a^2 }## term, or even both? More specifically, consider the following.

The integral is dominated by the value of the integrand at ##x = a## so focusing the terms in the product in a small interval ##\delta a ## around ##x = a##, we are basically dealing with a product of the form ##e^{-36 y^2}e^{252 y^3}## where ##y \equiv \frac{x - a}{a}\in [-\delta a/a, \delta a/a]##.
I've attached a plot of what each individual exponential looks like. It's clear then that in an extremely small neighborhood of the origin, which is where ##y## lies, the cubic exponential is almost entirely equal to unity so we can expand ##e^{252 y^3} \approx 1 + 252 y^3## to an excellent approximation.

The quadratic exponential, while not as close to unity in the origin neighborhood as the cubic term, still drops off very, very slowly from unity. In principle I suppose we could expand this to first order as well. Do we choose not to do it simply because it is an even better approximation to expand the cubic term at which point the quadratic term need not be expanded since the resulting integral can be evaluated analytically?

Don't know why the saddle point method gives a different result.

Perhaps the integrals are very sensitive to the method of approximation?

Thank you very much for the help by the way!

 

[TSny]

Well to second order we have [tex]\langle x \rangle = \frac{\int_{0}^{\infty} x e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} dx}{\int_{0}^{\infty} e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} dx} \approx \frac{\int_{a}^{\infty} x e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} dx}{\int_{a}^{\infty} e^{-36\beta U_0 \frac{(x - a)^2}{a^2}} dx}[/tex]

I don't understand your first approximation where you replace the lower limit of integration by ##a##.

For ##\beta U_0 \gg 1##, the exponential should damp out very quickly as you move away from x = ##a## to the left or right. So, to a very good approximation you can replace the lower limit by -∞. Then the two integrals can be done. I find ##<x> = a## independent of ##T##.

I understand why we can extend the integral to ##-\infty## to an excellent approximation but why is it valid to expand ##e^{252 \beta U_0 (x-a)^3/a^3 }## to first order? This is what I was originally stuck on actually. In the integral itself we have the product ##e^{252 \beta U_0 (x-a)^3/a^3 }e^{-36 \beta U_0 (x-a)^2/a^2 }##. So why expand the ## e^{252 \beta U_0 (x-a)^3/a^3 }## term and not the ##e^{-36 \beta U_0 (x-a)^2/a^2 }## term, or even both? More specifically, consider the following.

The integral is dominated by the value of the integrand at ##x = a## so focusing the terms in the product in a small interval ##\delta a ## around ##x = a##, we are basically dealing with a product of the form ##e^{-36 y^2}e^{252 y^3}## where ##y \equiv \frac{x - a}{a}\in [-\delta a/a, \delta a/a]##.
I've attached a plot of what each individual exponential looks like. It's clear then that in an extremely small neighborhood of the origin, which is where ##y## lies, the cubic exponential is almost entirely equal to unity so we can expand ##e^{252 y^3} \approx 1 + 252 y^3## to an excellent approximation.

The quadratic exponential, while not as close to unity in the origin neighborhood as the cubic term, still drops off very, very slowly from unity. In principle I suppose we could expand this to first order as well. Do we choose not to do it simply because it is an even better approximation to expand the cubic term at which point the quadratic term need not be expanded since the resulting integral can be evaluated analytically?

Perhaps the integrals are very sensitive to the method of approximation?

By approximating the exponential of the cubic part, you are left with integrals that can be done. I found it useful to plot ##e^{-\beta U(x)}## for various choices of U(x) and for some fixed choice of ##\beta U_0##, say 10:

U(x) = the exact Lennard-Jones potential
U(x) = the quadratic approximation = ##U_0 \left(-1+36(x -1)^2 \right)##
U(x) = the cubic approximation = ##U_0\left(-1+36(x -1)^2 -252(x-1)^3 \right)##
U(x) = the cubic approximation followed by letting ##e^{\;\beta U_0 252(x-1)^3} \approx 1+\beta U_0 252(x-1)^3##

Here, I've assumed measuring ##x## in units of ##a##.

It is interesting to see how the last approximation listed above gives a good fit to the exact ##e^{-\beta U(x)}## even when x is far from 1. However, the next to last approximation above causes ##e^{-\beta U(x)}## to "blow up" as you move too far away from x = 1 in the positive x direction.

 

[WannabeNewton]

Hi TSny! Sorry for the extremely late reply, I actually got stuck on a totally different problem haha (which you will probably see a thread on soon enough!)

I don't understand your first approximation where you replace the lower limit of integration by ##a##.

##\int_0^{\infty} \approx \int _a ^{\infty}##, is it not, since ##\int_0^a \approx 0##? Why is this not valid?

For ##\beta U_0 \gg 1##, the exponential should damp out very quickly as you move away from x = ##a## to the left or right. So, to a very good approximation you can replace the lower limit by -∞. Then the two integrals can be done. I find ##<x> = a## independent of ##T##.

I completely agree with you there and I think what you said agrees with what I did in the OP. But if that holds then certainly ##\int_0^{\infty} \approx \int _a ^{\infty} ## should hold as well since to the left of ##x = a## the integral is damped out very rapidly.

By approximating the exponential of the cubic part, you are left with integrals that can be done. I found it useful to plot ##e^{-\beta U(x)}## for various choices of U(x) and for some fixed choice of ##\beta U_0##, say 10:

U(x) = the exact Lennard-Jones potential
U(x) = the quadratic approximation = ##U_0 \left(-1+36(x -1)^2 \right)##
U(x) = the cubic approximation = ##U_0\left(-1+36(x -1)^2 -252(x-1)^3 \right)##
U(x) = the cubic approximation followed by letting ##e^{\;\beta U_0 252(x-1)^3} \approx 1+\beta U_0 252(x-1)^3##

Here, I've assumed measuring ##x## in units of ##a##.

It is interesting to see how the last approximation listed above gives a good fit to the exact ##e^{-\beta U(x)}## even when x is far from 1. However, the next to last approximation above causes ##e^{-\beta U(x)}## to "blow up" as you move too far away from x = 1 in the positive x direction.

Thank you. This was quite instructive. However what I still do not necessarily get is the motivation behind expanding the cubic as opposed to the quadratic. Is it simply because of what I said in the previous post i.e. that near ##x = a##, which is where the integral dominates, the cubic is as close to remaining at constant unity as possible without actually being unity whereas the quadratic drops off from unity albeit doing so extremely slowly?

 

[TSny]

This is just my opinion on it. So take it with a grain of salt.

One reason not to expand the quadratic exponential ##e^{-c*(x-a)^2}## is that the quadratic exponential makes the integrand get small for x not close to x = a. This allows the lower limit of integration to be changed from 0 to -∞. If you expanded the exponential of the quadratic expression about x = a and kept only 2 or 3 terms, this would not be true. You want an approximation that is good for all x in the phase space (not just very near x = a).

Another reason not to expand the quadratic exponential is that the integration of this type of exponential is a Gaussian type integral which can be easily carried out if integrating over all x. When the cubic exponential is expanded, you then get a polynomial times a Gaussian which is still not hard to integrate.

However, I'm with you in not being very comfortable with it! I cannot tell in general when this type of approximation magic is going to work well. It was only after graphing the various approximations that I could see that expanding the cubic exponential while not expanding the quadratic exponential appeared to be a good approximation. If the form of the potential was changed from the Lennard-Jones to some other form, I would need to graph the approximations again in order to have confidence in the procedure.

Sorry I can't give a good argument for why these tricks should work in a general case.

 

[WannabeNewton]

Thank you very much TSny. I think I understand for the most part. Plotting the exponential is definitely a very effective method for guessing an approximation so I'm with you on that.

However I'm still confused on this point:

##\int_0^{\infty} \approx \int _a ^{\infty}##, is it not, since ##\int_0^a \approx 0##? Why is this not valid?

I ask only because doing this gives that hideous expression for ##\alpha## to second order in post #7 whereas extending the lower bound to ##-\infty## just gives ##\alpha = 0## to second order.

 

[TSny]

Here's the graph of the integrand ##x e^{-\beta U}## for ##\beta U_0 = 5## using the exact Lennard-Jones potential. (Again, I chose units so ##a = 1##).

As you can see, you don't want to leave out the contribution from about ##.9 a## to ##a##.

 

[WannabeNewton]

Oh I see what you mean, thanks. On the other hand, with the extension to ##-\infty## we aren't leaving anything out, we're just including completely negligible terms.

But in light of that, are my arguments in the first half of post #1 for why we can actually Taylor expand the integrand in the first place about ##x = a## while still keeping the bounds ##[0,\infty)## (or ##(-\infty,\infty)## since the extension includes only negligible terms) actually valid?

 

[TSny]

If you look at the derivation of the saddle point result given in the link I previously posted, you can see that the derivation approaches the problem by keeping the quadratic exponential but expanding the higher order exponential. But note that the expansion is in terms of 1/√A, where A would correspond to U_o/(kT) in the Lennard Jones problem.

I decided to check the expression for C₂ given below equation (1.7). I do not get the same numerical coefficients in the expression shown there. Instead, I get $$C_2 = -\frac{f''}{2fg''}+\frac{f'g'''}{2f(g'')^2} +\frac{g''''}{8(g'')^2} -\frac{5(g''')^2}{24(g'')^3} $$

When I use this expression to find the thermal expansion coefficient, I get $$\alpha \approx \frac{7}{48} \frac{k_B}{U_0} \approx 0.15 \frac{k_B}{U_0} $$

 

[WannabeNewton]

Sorry I'm not following. What did you do differently in your derivation of ##C_2## from that of the saddle-point approximation paper?

 

[TSny]

I followed the method given in the link.

If you look at equation (1.7) you see that the coefficients ##C_{2n}## come from carrying out the integrations in the first line of (1.7) and collecting like terms in ##\large \frac{1}{A^n}##.

When I did the integration I found different numerical coefficients in the resultant terms for ##C_2##.

If you take the result for ##C_2## as given in the link, and everywhere replace ##g''(x_0)## by ##g''(x_0)/2##, then you get the result that I got.

So, I'm thinking that the author of the link might have left out the 1/2 factor in the exponential ##\large e^{y^2g''(x_0)/2}## in the first line of (1.7) when doing the integration.

Otherwise, I'm the one making an error. But it's interesting that I now get the 0.15 result.

 

[WannabeNewton]

I did the integration by hand as well and got exactly what you got so perhaps the author did make a typo. And indeed the corrected coefficient gives the desired value of ##.15##. Just to be safe I'll do the calculation following the discussion in Arfken and Weber for the saddle-point approximation and see if I get the ##.15## again. Thanks so much for the help!