Advanced Chemistry - Mono Mole

April 30, 2024May 8, 2025

Taylor series of single-variable and multi-variable functions

The Taylor series is a way to express a function as an infinite sum of terms, each of which is derived from the function’s derivative value at a specific point $a$ .

Consider the power series $f(x)=c_0+c_1(x-a)+c_2(x-a)^2+\cdots+c_n(x-a)^n$ . We have $f(a)=c_0$ , $f'(a)=c_1$ , $f''(a)=2c_2$ , $f'''(a)=6c_3$ and so on. Therefore,

$f(x)=f(a)+\frac{f'(a)}{1!}(x-a)+\frac{f''(a)}{2!}(x-a)^2+\cdots+\frac{f^n(a)}{n!}(x-a)^n\;\;\;\;\;\;\;\;32b$

Eq32b is called a Taylor series, which approximates a function in the neighbourhood of the point $a$ . When $a=0$ , the Taylor series is also known as a Maclaurin series:

$f(x)=f(0)+\frac{f'(0)}{1!}x+\frac{f''(0)}{2!}x^2+\cdots+\frac{f^n(0)}{n!}x^n$

The input value $x$ of $f$ in eq32b represents points in the domain of $f$ that are near $a$ . In other words, we can express $f(x)$ as $f(a+t(x-a))$ , where $a$ is a constant, $t$ is a variable, and $x-a$ represents a small change in $a$ that is scaled by $t$ . This alternate expression of $f$ is useful when dealing with multiple-variable functions.

For a multivariable function $f(x_1,x_2,\cdots,x_n)$ , the Taylor series about the point $(a_1,a_2,\cdots,a_n)$ is

$\begin{align}f(x_1,x_2,\cdots,x_n)&=f(a_1,a_2,\cdots,a_n)+\sum_{i=1}^n\frac{\partial f(a_1,a_2,\cdots,a_n)}{\partial x_i}(x_i-a_i)\\&+\frac{1}{2!}\frac{\partial^2 f(a_1,a_2,\cdots,a_n)}{\partial x_i\partial x_j}(x_i-a_i)(x_j-a_j)+\cdots\;\;\;\;\;\;\;\;32c\end{align}$

Question

How do we use vector notation to express the function $f(x_1,x_2,\cdots,x_n)=c_1x_1+c_2x_2+\cdots+c_nx_n$ ?

Answer

If we let $\hat{x}=(x_1,x_2,\cdots,x_n)$ , then we write $f(\hat{x})$ in place of $f(x_1,x_2,\cdots,x_n)$ . In other words, we have $f(\hat{x})=\hat{c}\cdot\hat{x}$ , where $\hat{c}=(c_1,c_2,\cdots,c_n)$ . Consequently, $f$ can be regarded as

1. a function of $n$ variables $x_1,x_2,\cdots,x_n$ ,
2. a function of a single point variable $(x_1,x_2,\cdots,x_n)$ , or
3. a function of a single vector variable $\hat{x}$ .

To derive eq32c, let $\hat{x}=(x_1,x_2,\cdots,x_n)$ and $\hat{a}=(a_1,a_2,\cdots,a_n)$ . Consider the function

$g(t)=f(\hat{b}=\hat{a}+t(\hat{x}-\hat{a}))\;\;\;\;\;\;\;\;32d$

which implies that $\hat{b}=(x_1,x_2,\cdots,x_n)$ .

$f(\hat{a}+t(\hat{x}-\hat{a}))$ is a multivariable function, whose input is the vector $\hat{a}+t(\hat{x}-\hat{a})$ , which varies with $t$ . In the domain of $f(x_1,x_2,\cdots,x_n)$ , $\hat{a}$ is the vector representing a point where $f$ is expanded, and $(\hat{x}-\hat{a})$ is a fixed vector that determines the direction of the displacement from the point $\hat{a}$ . In other words, $(\hat{x}-\hat{a})$ is not a variable of $g(t)$ , but a parameter that we choose before plotting $g(t)$ . The function $g(t)$ is therefore a single variable function of $t$ , and its points are of the form $(t,g(t))$ .

The Maclaurin series expansion of $g(t)$ when $t=1$ is

$g(1)=g(0)+\frac{g'(0)}{1!}+\frac{g''(0)}{2!}+\cdots\;\;\;\;\;\;\;\;32e$

According to the multivariable chain rule, ${\frac{df(\hat{b})}{dt}=\frac{\partial f(\hat{b})}{\partial x_1}\frac{dx_1}{dt}+\frac{\partial f(\hat{b})}{\partial x_2}\frac{dx_2}{dt}+\cdots+\frac{\partial f(\hat{b})}{\partial x_n}\frac{dx_n}{dt}=\nabla f(\hat{b})\cdot\frac{d\hat{x}}{dt}}$ , where ${\nabla=\boldsymbol{\mathit{i}}\frac{\partial}{\partial x_1}+\boldsymbol{\mathit{j}}\frac{\partial}{\partial x_2}+\cdots+\boldsymbol{\mathit{k}}\frac{\partial}{\partial x_n}}$ and $\hat{x}=\boldsymbol{\mathit{i}}x_1+\boldsymbol{\mathit{j}}x_2+\cdots+\boldsymbol{\mathit{k}}x_n$ . So,

$g'(t)=\nabla f(\hat{b})\cdot\frac{d\hat{x}}{dt}=\cdots=\sum_{i=1}^n\frac{\partial f(\hat{b})}{\partial x_i}(x_i-a_i)\;\;\;\;\;\;\;\;32f$

The second derivative $g''(t)=f''(\hat{b})$ is

$g''(t)=\frac{d^2f(\hat{b})}{dt^2}=\frac{d}{dt}\biggr\[\frac{\partial f(\hat{b})}{\partial x_1}\frac{dx_1}{dt}+\frac{\partial f(\hat{b})}{\partial x_2}\frac{dx_2}{dt}+\cdots+\frac{\partial f(\hat{b})}{\partial x_n}\frac{dx_n}{dt}\biggr\]$

Using the multivariable chain rule again, we get

$g''(t)=\frac{d^2f(\hat{b})}{dt^2}=\sum_{i,j=1}^n\frac{\partial^2 f(\hat{b})}{\partial x_i\partial x_j}\frac{dx_i}{dt}\frac{dx_j}{dt}+\sum_{k=1}^n\frac{\partial f(\hat{b})}{\partial x_k}\frac{d^2x_k}{dt^2}\;\;\;\;\;\;\;\;32g$

With reference to eq32d, $\hat{b}=\hat{a}+t(\hat{x}-\hat{a})$ with components $x_i=a_i+t(x_i-a_i)$ . It follows that $\frac{dx_i}{dt}=x_i-a_i$ , $\frac{dx_j}{dt}=x_j-a_j$ and $\frac{d^2x_k}{dt^2}=0$ . Eq32g becomes

$g''(t)=\frac{d^2f(\hat{b})}{dt^2}=\sum_{i,j=1}^n\frac{\partial^2 f(\hat{b})}{\partial x_i\partial x_j}(x_i-a_i)(x_j-a_j)\;\;\;\;\;\;\;\;32h$

Furthermore, when $t=0$ , we have $\hat{b}=\hat{a}$ . So, eq32f and eq32h become

$g'(0)=\sum_{i=1}^n\frac{\partial f(\hat{a})}{\partial x_i}(x_i-a_i)\;\;\;\;\;\;\;\;32i$

and

$g''(0)=\sum_{i,j=1}^n\frac{\partial^2 f(\hat{a})}{\partial x_i\partial x_j}(x_i-a_i)(x_j-a_j)\;\;\;\;\;\;\;\;32j$

respectively.

Substituting eq32i and eq32j in eq32e and noting that $g(0)=f(\hat{a})$ and $g(1)=f(\hat{x})$ , we have eq32c. Finally, the Maclaurin series of a multivariable function is

$\begin{align}f(x_1,x_2,\cdots,x_n)&=f(a_1=0,a_2=0,\cdots,a_n=0)\\&+\sum_{i=1}^n\frac{\partial f(a_1=0,a_2=0,\cdots,a_n=0)}{\partial x_i}x_i\\&+\frac{1}{2!}\frac{\partial^2 f(a_1=0,a_2=0,\cdots,a_n=0)}{\partial x_i\partial x_j}x_ix_j\\&+\cdots\;\;\;\;\;\;\;\;32k\end{align}$

Next article: Generating function and Rodrigues’ formula for the Hermite polynomials

Previous article: Recurrence relations of the Hermite polynomials

Content page of vibrational spectroscopy

Content page of advanced chemistry

Main content page

April 30, 2024August 30, 2024

The quantum harmonic oscillator

The quantum harmonic oscillator is a model for studying an atom that moves back and forth about an equilibrium point. The model’s central premise involves the use of the potential energy term of the classical harmonic oscillator to construct the Schrodinger equation for the harmonic oscillator:

$\biggr$\frac{p^2}{2m}+\frac{1}{2}kx^2\biggr$\psi(x)=E\psi(x)\;\;\;\;\;\;\;\;4$
which is equivalent to

$\biggr$-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+\frac{1}{2}kx^2\biggr$\psi(x)=E\psi(x)$

where $x$ is the displacement of the mass from its equilibrium position.

Question

Show that $\frac{p^2}{2m}$ is equivalent to $-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}$ .

Answer

The de Broglie relation is $p=\hbar k$ , where $k=\frac{2\pi}{\lambda}$ . The derivative of the wavefunction $\psi=Ae^{ikx}$ with respect to $x$ is

$\frac{d\psi}{dx}=ikAe^{ikx}=ik\psi$

$\frac{\hbar}{\sqrt{2m}}\frac{d}{dx}\psi=i\frac{\hbar}{\sqrt{2m}}k\psi=i\frac{p}{\sqrt{2m}}\psi$

This suggests that the operator $\frac{\hbar}{\sqrt{2m}}\frac{d}{dx}$ is equivalent to $i\frac{p}{\sqrt{2m}}$ , and hence, $\hat{p}=\frac{\hbar}{i}\frac{d}{dx}$ . Squaring the momentum operator $\hat{p}$ and dividing by $2m$ gives $-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}$ .

Substituting $k=m\omega^2$ in the above equation and rearranging, we have

$\frac{d^2}{dx^2}\psi(x)+\biggr$\frac{2m}{\hbar^2}E-\frac{m^2\omega^2}{\hbar^2}\biggr$\psi(x)=0\;\;\;\;\;\;\;\;5$

The probability of locating an electron in the Hilbert space must be finite, i.e. $\int_{-\infty}^{\infty}\vert\psi(x)\vert^2dx<\infty$ . This implies that $\psi(x)$ vanishes as $x\rightarrow\pm\infty$ . We call such a well-behaved wavefunction square-integrable. Hence, we need to study the asymptotic characteristics of eq5 to find possible solutions.

Let’s begin by simplifying the differential equation using a change of variable. Substituting $x=y\biggr$\frac{\hbar}{m\omega}\biggr$^{1/2}$ and $\epsilon=\frac{E}{\hbar\omega}$ in eq5, we have

$\frac{d^2\psi(y)}{dy^2}+\(2\epsilon-y^2)\psi(y)=0\;\;\;\;\;\;\;\;6$

Since $2\epsilon$ is a constant, $y^2\gg 2\epsilon$ as $y\rightarrow\pm\infty$ , and eq6 approximates to

$\frac{d^2\psi(y)}{dy^2}-y^2\psi(y)=0\;\;\;\;\;\;\;\;7$

This suggests that the solution to eq6 could be a Gaussian function because $\psi(y)=e^{-\frac{y^2}{2}}$ is a solution to eq7. Let’s try $\psi(y)=u(y)e^{-\frac{y^2}{2}}$ , where $u(y)$ is a function of $y$ . Substituting $\psi(y)=u(y)e^{-\frac{y^2}{2}}$ in eq6 and computing the derivatives, we have

$\frac{d^2u(y)}{dy^2}-2y\frac{du(y)}{dy}+(2\epsilon-1)u(y)=0\;\;\;\;\;\;\;\;8$

Eq8 is known as the Hermite differential equation.

Question

Show that we can use a power series to solve eq8 around $y=0$ .

Answer

When $y=0$ , eq8 simplifies to $\frac{d^2u(y)}{dy^2}+c^2u(y)=0$ , where $c=\sqrt{2\epsilon-1}$ . It can be easily verified that $u=Acos(cy)$ and $u=Bsin(cy)$ are solutions to $\frac{d^2u(y)}{dy^2}+c^2u(y)=0$ . As the Taylor series for $u=Acos(cy)$ and $u=Bsin(cy)$ are $u=A\sum_{k=0}^{\infty}(-1)^k\frac{c^{2k}x^{2k}}{(2k)!}$ and $u=B\sum_{k=0}^{\infty}(-1)^k\frac{c^{2k+1}x^{2k+1}}{(2k+1)!}$ , respectively, a power series is a solution to eq8 around $y=0$ .

The above Q&A suggests that a possible solution to eq8 is a power series for a range of values of $y$ . Let $u(y)=\sum_{n=0}^{\infty}C_ny^n$ and substitute $\frac{du(y)}{dy}=\sum_{n=0}^{\infty}nC_ny^{n-1}$ and $\frac{d^2u(y)}{dy^2}=\sum_{n=0}^{\infty}n(n-1)C_ny^{n-2}$ in eq8 to give

$\sum_{n=0}^{\infty}n(n-1)C_ny^{n-2}-2y\sum_{n=0}^{\infty}nC_ny^{n-1}+(2\epsilon-1)\sum_{n=0}^{\infty}C_ny^{n}=0$

Substituting $2y\sum_{n=0}^{\infty}nC_ny^{n-1}=2\sum_{n=0}^{\infty}nC_ny^n$ and $\sum_{n=0}^{\infty}n(n-1)C_ny^{n-2}=\sum_{n=0}^{\infty}(n+2)(n+1)C_{n+2}y^n$ in the above equation, we have

$\sum_{n=0}^{\infty}\[(n+2)(n+1)C_{n+2}-2nC_n+(2\epsilon-1)C_n\]y^n=0$

The above equation is only true for all values of $y$ if every coefficient for each power of $y$ is zero. So, $(n+2)(n+1)C_{n+2}-2nC_n+(2\epsilon-1)C_n=0$ , which rearranges to

$C_{n+2}=\frac{2n+1-2\epsilon}{(n+2)(n+1)}C_n\;\;\;\;\;\;\;\;9$

Eq9 is a recurrence relation. If we know the value of $C_0$ , we can use the relation to find $C_2,C_4,C_6,\cdots$ . Similarly, if we know $C_1$ , we can find $C_3,C_5,C_7,\cdots$ .

Comparing the recurrence relations for even-labelled coefficients,

$C_{2n}=-\[2\epsilon-1-2(2n-2)\]\frac{(2n-2)!}{(2n)!}C_{2n-2}\;\;\;where\;n=1,2,3\cdots\;\;\;\;10$

Similarly, the recurrence relations for odd-labelled coefficients can be expressed as

$C_{2n+1}=-\[2\epsilon-3-2(2n-2)\]\frac{(2n-1)!}{(2n+1)!}C_{2n-1}\;\;\;where\;n=1,2,3\cdots\;\;\;\;11$

Rewriting $\psi(y)$ as a sum of two power series, one containing terms with even-labelled coefficients and the other with odd-labelled coefficients,

$\psi(y)=\biggr\{C_0+\sum_{l=1}^{\infty}\biggr\[-\[2\epsilon-1-2(2l-2)\]\frac{C_{2l-2}}{\frac{(2l)!}{(2l-2)!}}\biggr\]y^{2l}\biggr\}e^{-\frac{y^2}{2}}+\\\biggr\{C_1y+\sum_{l=1}^{\infty}\biggr\[-\[2\epsilon-3-2(2l-2)\]\frac{C_{2l-1}}{\frac{(2l+1)!}{(2l-1)!}}\biggr\]y^{2l+1}\biggr\}e^{-\frac{y^2}{2}}\;\;\;\;\;\;\;\;12$

According to L’Hopital’s rule, $\psi(y)$ appears to be a possible solution to eq7 because each of the terms in eq12 (simplified to $Cy^n/e^{\frac{y^n}{2}}$ ) approaches zero for large $y$ :

$\lim\limits_{y\rightarrow\pm\infty}\frac{Cy^n}{e^{\frac{y^2}{2}}}=\lim\limits_{y\rightarrow\pm\infty}\frac{nCy^{n-1}}{e^{\frac{y^2}{2}}}=\cdots=\lim\limits_{y\rightarrow\pm\infty}\frac{k}{e^{\frac{y^2}{2}}}=0$

where $k$ is a constant.

However, this does not guarantee that the series as a whole also converges because the sum of an infinite number of infinitesimal terms can still be infinite, as is the case for a harmonic series. To see how the two infinite series behave for large $y$ , we carry out the ratio test as follows:

where we have set $n=2l$ and $n=2l+1$ in eq9 for the even series and odd series respectively.

By comparison, the ratio test for the power series expansion of $e^{y^2}=\sum_{l=0}^{\infty}\frac{y^{2l}}{l!}$ is also $\frac{1}{l}$ . Therefore, both infinite series behave the same as $e^{y^2}$ and diverge for large $y$ . To ensure that $\psi(y)$ is square-integrable, we need to truncate either one of the series after some finite terms $l$ and let all the coefficients of the other series be zero.

Question

Can we instead truncate both series of eq12 after some $l$ ?

Answer

Since the value of $l$ is arbitrary, the sum of two truncated series, each with finite terms, may still have an infinite number of terms. The only way to guarantee that $u(y)$ has finite terms is to truncate one of the series and let $C_0$ , if the odd series is truncated, or $C_1$ , if the even series is truncated, be zero.

To truncate either series, we let the numerator of eq9 be zero so that every successive term in the selected series is zero as well. This implies that

$\epsilon=n+\frac{1}{2}\;\;\;\;\;\;\;\;13$

Let’s rewrite eq12 as

$\psi(y)=\biggr\{\begin{matrix}\biggr\{C_0+\sum_{l=1}^{m}\biggr\[-\[2\epsilon-1-2(2l-2)\]\frac{C_{2l-2}}{\frac{(2l)!}{(2l-2)!}}\biggr\]y^{2l}\biggr\}e^{-\frac{y^2}{2}}&even\;series\;l=1,2,\cdots\\\biggr\{C_1y+\sum_{l=1}^{m}\biggr\[-\[2\epsilon-3-2(2l-2)\]\frac{C_{2l-1}}{\frac{(2l+1)!}{(2l-1)!}}\biggr\]y^{2l+1}\biggr\}e^{-\frac{y^2}{2}}&odd\;series\;l=1,2,\cdots\end{matrix}\;\;\;\;14$

Since eq14 is the result of substituting $\psi(y)=u(y)e^{-\frac{y^2}{2}}$ in eq6, it is a solution to eq6.

Question

Show that $\psi_4(y)=C_0\biggr\[1-\frac{2n}{2!}y^2+\frac{2^2n(n-2)}{4!}y^4\biggr\]e^{-\frac{y^2}{2}}$ .

Answer

We expand the even series in eq14 as follows:

$\psi_4(y)=\biggr\[C_0-(2\epsilon-1)\frac{C_0}{2!}y^2-2(2\epsilon-1-4)\frac{C_2}{4!}y^4\biggr\]e^{-\frac{y^2}{2}}$

Substituting eq9, where $C_2=\frac{1-2\epsilon}{2}C_0$ , in the above equation, we have

$\psi_4(y)=\biggr\[C_0-(2\epsilon-1)\frac{C_0}{2!}y^2-\frac{(2\epsilon-1-4)(1-2\epsilon)}{4!}C_0y^4\biggr\]e^{-\frac{y^2}{2}}$

Substituting eq13 in the above equation gives the required expression.

Re-expanding the even series of eq14 using eq9 and eq13, we have

$\boldsymbol{\mathit{n}}$	$\boldsymbol{\mathit{\psi_n(y)}}$
$0$	$C_0e^{-\frac{y^2}{2}}$
$2$	$C_0\biggr\[1-\frac{2n}{2!}y^2\biggr\]e^{-\frac{y^2}{2}}$
$4$	$C_0\biggr\[1-\frac{2n}{2!}y^2+\frac{2^2n(n-2)}{4!}y^4\biggr\]e^{-\frac{y^2}{2}}$
$6$	$C_0\biggr\[1-\frac{2n}{2!}y^2+\frac{2^2n(n-2)}{4!}y^4-\frac{2^3n(n-2)(n-4)}{6!}y^6\biggr\]e^{-\frac{y^2}{2}}$

Comparing the coefficients of each even series in the above table, we have

$\psi_n(y)=C_0\biggr\[1-\frac{2n}{2!}y^2+\frac{2^2n(n-2)}{4!}y^4+\cdots+\frac{2^n$\frac{n}{2}$!}{(-1)^{\frac{n}{2}}n!}y^n\biggr\]e^{-\frac{y^2}{2}}\;\;\;\;\;\;\;\;15a$

where $n=0,2,4,\cdots$ .

Similarly, the re-expansion of the odd series of eq14 using eq9 and eq13 gives

$\boldsymbol{\mathit{n}}$	$\boldsymbol{\mathit{\psi_n(y)}}$
$1$	$C_1ye^{-\frac{y^2}{2}}$
$3$	$C_1\biggr\[y-\frac{2(n-1)}{3!}y^3\biggr\]e^{-\frac{y^2}{2}}$
$5$	$C_1\biggr\[y-\frac{2(n-1)}{3!}y^3+\frac{2^2(n-1)(n-3)}{5!}y^5\biggr\]e^{-\frac{y^2}{2}}$
$7$	$C_1\biggr\[y-\frac{2(n-1)}{3!}y^3+\frac{2^2(n-1)(n-3)}{5!}y^5-\frac{2^3(n-1)(n-3)(n-5)}{7!}y^7\biggr\]e^{-\frac{y^2}{2}}$

Comparing the coefficients of each odd series in the above table, we have

$\psi_n(y)=C_1\biggr\[y-\frac{2(n-1)}{3!}y^3+\frac{2^2(n-1)(n-3)}{5!}y^5+\cdots+\frac{2^n$\frac{n-1}{2}$!}{(-1)^{\frac{n-1}{2}}2(n!)}y^n\biggr\]e^{-\frac{y^2}{2}}\;\;\;\;\;\;\;\;15b$

where $n=1,3,5,\cdots$ .

The recurrence relation of eq9 that led to eq14 only defines the value of $C_n$ in terms of $C_{n+2}$ . This implies that there are many solution sets depending on the value of $C_{n+2}$ . The convention is to set the leading coefficients (i.e. the coefficients of $y^n$ ) of eq15a and eq15b, which are independent of each other, as $2^n$ . Therefore,

$C_0\frac{2^n$\frac{n}{2}$!}{(-1)^{\frac{n}{2}}n!}=2^n\;\;or\;\;C_0=(-1)^{\frac{n}{2}}\frac{n!}{(\frac{n}{2})!}\;\;\;\;\;\;\;\;16$

and

$C_1\frac{2^n$\frac{n-1}{2}$!}{(-1)^{\frac{n-1}{2}}2(n!)}=2^n\;\;or\;\;C_1=(-1)^{\frac{n-1}{2}}\frac{2(n!)}{(\frac{n-1}{2})!}\;\;\;\;\;\;\;\;17$

Substituting eq16 and eq17 in eq14, we have $\psi(y)=H_n(y)e^{-\frac{y^2}{2}}$ , where $H_n(y)$ are the Hermite polynomials:

$H_n(y)=\biggr\{\begin{matrix}\biggr\{(-1)^{\frac{n}{2}}\frac{n!}{(\frac{n}{2})!}+\sum_{l=1}^{m}\biggr\[-\[2\epsilon-1-2(2l-2)\]\frac{C_{2l-2}}{\frac{(2l)!}{(2l-2)!}}\biggr\]y^{2l}\biggr\}e^{-\frac{y^2}{2}}&even\;series\;l=1,2,\cdots\\\biggr\{(-1)^{\frac{n-1}{2}}\frac{2(n!)}{(\frac{n-1}{2})!}y+\sum_{l=1}^{m}\biggr\[-\[2\epsilon-3-2(2l-2)\]\frac{C_{2l-1}}{\frac{(2l+1)!}{(2l-1)!}}\biggr\]y^{2l+1}\biggr\}e^{-\frac{y^2}{2}}&odd\;series\;l=1,2,\cdots\end{matrix}\;\;\;\;18$

The first few Hermite polynomials, un-normalised eigenfunctions of eq6 and their corresponding eigenvalues are

When $y=0$ , we can express the Hermite polynomials for odd $n$ as

$H_{2n+1}(0)=0\;\;\;n=0,1,2,\cdots\;\;\;\;\;\;\;\;19$

Similarly, when $y=0$ , we can express the Hermite polynomials for even $n$ as

$H_{2n}(0)=\frac{(-1)^n(2n)!}{n!}\;\;\;n=0,1,2,\cdots\;\;\;\;\;\;\;\;20$

Eq19 and eq20 are useful in proving the generating function for Hermite polynomials.

The table above also shows that the energy states $E_n$ of a quantum-mechanical harmonic oscillator are quantised:

$E_n=\epsilon\hbar\omega=\biggr$n+\frac{1}{2}\biggr$\hbar\omega\;\;\;n=0,1,2,\cdots\;\;\;\;\;21$

with the ground state $E_0=\frac{1}{2}\hbar\omega$ being non-zero. This energy is called the zero-point energy of the harmonic oscillator.

Question

Using the uncertainty principle, explain why the ground state of a quantum harmonic oscillator is non-zero.

Answer

If $E_0=0$ , then both the kinetic energy $E_k=\frac{p^2}{2m}$ and the potential energy $V$ of the oscillator are zero. $E_k=0$ implies that the momentum $p$ is zero and $\Delta p=0$ . $V=0$ means that the atom is at the equilibrium position, where $\Delta x=0$ . This would violate the uncertainty principle.

Since $H_n(y)$ is a solution to eq8, we can rewrite eq8 as

$\frac{d^2H_n(y)}{dy^2}-2y\frac{dH_n(y)}{dy}+(2\epsilon-1)H_n(y)=0\;\;\;\;\;\;\;\;22$

To complete the prove that $\psi_n(y)=H_n(y)e^{-\frac{y^2}{2}}$ is a solution to eq6, we will derive the normalisation constant $N_n$ of $\psi_n(y)=N_nH_n(y)e^{-\frac{y^2}{2}}$ and prove that the wavefunctions are orthogonal to one another in the next few articles.

Question

Do we have to change the variable $y$ back to $x$ ?

Answer

The Hermite polynomials are usually expressed in terms of $y$ because $H_n\biggr$\sqrt{\frac{m\omega}{\hbar}}x\biggr$$ looks more complicated:

Hermite polynomials

$H_0\biggr$\sqrt{\frac{m\omega}{\hbar}}x\biggr$=1$

$H_1\biggr$\sqrt{\frac{m\omega}{\hbar}}x\biggr$=2\sqrt{\frac{m\omega}{\hbar}}x$

$H_2\biggr$\sqrt{\frac{m\omega}{\hbar}}x\biggr$=4\biggr$\sqrt{\frac{m\omega}{\hbar}}x\biggr$^2-2$

$\vdots$
However, the normalisation constant is defined as an integral with respect to $x$ because $x$ is the variable representing position in the one-dimensional space that the oscillator is moving in.

Next article: Orthogonality of the wavefunctions of the quantum harmonic oscillator

Previous article: The harmonic oscillator

Content page of vibrational spectroscopy

Content page of advanced chemistry

Main content page

April 30, 2024

Orthogonality of the wavefunctions of the quantum harmonic oscillator

The orthogonality of the wavefunctions of the quantum harmonic oscillator $\psi_n(y)=H_ne^{-\frac{y^2}{2}}$ can be proven using the Hermite differential equation.

Substituting eq13 in eq22, we have

$\frac{d^2H_n}{dy^2}-2y\frac{dH_n}{dy}+2nH_n=0\;\;\;\;\;\;\;\;23$

where $H_n$ are the Hermite polynomials.

Question

Show that eq23 can be written as

$e^{y^2}\frac{d}{dy}\biggr$e^{-y^2}\frac{dH_n}{dy}\biggr$+2nH_n=0\;\;\;\;\;\;\;\;24$

Answer

Carrying out the derivative, the LHS of eq24 becomes

$e^{y^2}\biggr\[e^{-y^2}\frac{d^2H_n}{dy^2}+\frac{dH_n}{dy}(-2y)e^{-y^2}\biggr\]+2nH_n=\frac{d^2H_n}{dy^2}-2y\frac{dH_n}{dy}+2nH_n$

Multiplying eq24 (with a change of index from $n$ to $m$ ) by $H_n$ and subtracting the result from the product of $H_m$ and eq24, we have

${H_me^{y^2}\frac{d}{dy}\biggr$e^{-y^2}\frac{dH_n}{dy}\biggr$+2nH_mH_n-\biggr\[H_ne^{y^2}\frac{d}{dy}\biggr$e^{-y^2}\frac{dH_m}{dy}\biggr$+2mH_mH_n\biggr\]=0}$

$e^{-y^2}2(m-n)H_mH_n=H_m\frac{d}{dy}\biggr$e^{-y^2}\frac{dH_n}{dy}\biggr$-H_n\frac{d}{dy}\biggr$e^{-y^2}\frac{dH_m}{dy}\biggr$\;\;\;\;\;\;\;\;25$

Substituting ${H_m\frac{d}{dy}\biggr$e^{-y^2}\frac{dH_n}{dy}\biggr$=\frac{d}{dy}\biggr$H_me^{-y^2}\frac{dH_n}{dy}\biggr$-e^{-y^2}\frac{dH_n}{dy}\frac{dH_m}{dy}}$ and ${H_n\frac{d}{dy}\biggr$e^{-y^2}\frac{dH_m}{dy}\biggr$=\frac{d}{dy}\biggr$H_ne^{-y^2}\frac{dH_m}{dy}\biggr$-e^{-y^2}\frac{dH_m}{dy}\frac{dH_n}{dy}}$ in eq25 gives

$e^{-y^2}2(m-n)H_mH_n=\frac{d}{dy}\biggr\[e^{-y^2}\biggr$H_m\frac{dH_n}{dy}-H_n\frac{dH_m}{dy}\biggr$\biggr\]$

Integrating both sides of the above equation with respect to $y$ ,

$2(m-n)\int_{-\infty}^{\infty}e^{-y^2}H_mH_ndy=\int_{-\infty}^{\infty}\frac{d}{dy}\biggr\[e^{-y^2}\biggr$H_m\frac{dH_n}{dy}-H_n\frac{dH_m}{dy}\biggr$\biggr\]dy=0$

If $m\neq n$ , then $\int_{-\infty}^{\infty}e^{-y^2}H_mH_ndy=0$ . Since $\psi_n(y)=H_ne^{-\frac{y^2}{2}}=\psi^*_n(y)$ ,

$\int_{-\infty}^{\infty}\psi^*_m(y)\psi_n(y)dy=0\;\;\;\;m\neq n\;\;\;\;\;\;\;\;25a$

Next article: Recurrence relations of the Hermite polynomials

Previous article: The quantum harmonic oscillator

Content page of vibrational spectroscopy

Content page of advanced chemistry

Main content page

April 30, 2024August 4, 2025

Normalisation constant for the wavefunctions of the quantum harmonic oscillator

The normalisation constant ensures that the wavefunctions of the quantum harmonic oscillator are properly scaled, thereby maintaining the probabilistic interpretation of quantum states.

To derive the normalisation constant for the wavefunctions of the quantum harmonic oscillator, we begin with the multiplication of eq33 by $e^{-y^2}$ and by itself (with a change of the dummy variable $t$ to $s$ and the dummy index $n$ to $m$ ) to give

$e^{-y^2}e^{-t^2+2ty}e^{-s^2+2sy}=e^{-y^2}\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\frac{H_m(y)}{m!}\frac{H_n(y)}{n!}s^mt^n\;\;\;\;\;\;\;\;41$

Integrating with respect to $y$ ,

$\int_{-\infty}^{\infty}e^{-y^2-t^2+2ty-s^2+2sy}dy=\int_{-\infty}^{\infty}e^{-y^2}\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\frac{H_m(y)}{m!}\frac{H_n(y)}{n!}s^mt^ndy$

Substituting the orthogonal relation $\int_{-\infty}^{\infty}e^{-y^2}H_m(y)H_n(y)dy=0$ when $m\neq n$ in the above equation and simplifying, we have

$e^{2st}\int_{-\infty}^{\infty}e^{-(y-s-t)^2}dy=\sum_{n=0}^{\infty}\frac{(st)^n}{n!n!}\int_{-\infty}^{\infty}e^{-y^2}\[H_n(y)\]^2dy\;\;\;\;\;\;\;\;42$

Question

Show that $I=\int_{-\infty}^{\infty}e^{-a(y-s-t)^2}dy=\sqrt{\frac{\pi}{a}}$ , where $a$ is a constant.

Answer

Let $X=y-s-t$ and $dX=dy$ . So, $I=\int_{-\infty}^{\infty}e^{-aX^2}dX$ with the limits unchanged. At the same time, let $I=\int_{-\infty}^{\infty}e^{-aY^2}dY$ , where we have changed the dummy variable $X$ to $Y$ .

$I^2=\int_{-\infty}^{\infty}e^{-aX^2}dX\int_{-\infty}^{\infty}e^{-aY^2}dY=\int\int_{-\infty}^{\infty}e^{-a(X^2+Y^2)}dXdY$

In polar coordinates, $r^2=X^2+Y^2,X=rcos\theta,Y=rsin\theta$ . For infinitesimal changes, the change in area $dA$ can be approximated as a rectangle with sides $r$ and $rd\theta$ (see diagram below), with $dXdY=dA=rd\theta dr$ .

So,

$I^2=\int_{0}^{2\pi}d\theta\int_{0}^{\infty}re^{-ar^2}dr=2\pi\int_{0}^{\infty}re^{-ar^2}dr$

Let $u=r^2\;\Rightarrow\;du=2rdr$

$I^2=2\pi\int_{0}^{\infty}re^{-au}\frac{du}{2r}=\pi\int_{0}^{\infty}e^{-au}du=\frac{\pi}{a}$

$I=\int_{-\infty}^{\infty}e^{-a(y-s-t)^2}dy=\sqrt{\frac{\pi}{a}}\;\;\;\;\;\;\;\;43$

Using eq43 and the Maclaurin series expansion of $e^x=\sum_{n=0}^{\infty}\frac{x^n}{n!}$ , eq42 becomes

$\sqrt{\pi}\sum_{n=0}^{\infty}\frac{2^n(st)^n}{n!}=\sum_{n=0}^{\infty}\frac{(st)^n}{n!n!}\int_{-\infty}^{\infty}e^{-y^2}\[H_n(y)\]^2dy\;\;\;\;\;\;\;\;44$

Expanding eq44 and comparing the coefficients of $(st)^n$ ,

$n=0$	$\int_{-\infty}^{\infty}e^{-y^2}\[H_n(y)\]^2dy=\sqrt{\pi}$
$n=1$	$\int_{-\infty}^{\infty}e^{-y^2}\[H_n(y)\]^2dy=2\sqrt{\pi}$
$n=2$	$\int_{-\infty}^{\infty}e^{-y^2}\[H_n(y)\]^2dy=8\sqrt{\pi}$
$n=3$	$\int_{-\infty}^{\infty}e^{-y^2}\[H_n(y)\]^2dy=48\sqrt{\pi}$
$\vdots$	$\vdots$
$n=n$	$\int_{-\infty}^{\infty}e^{-y^2}\[H_n(y)\]^2dy=2^nn!\sqrt{\pi}$

In an earlier article, we have made a change of variable of ${x=y\biggr$\frac{\hbar}{m\omega}\biggr$^{\frac{1}{2}}}$ in deriving the un-normalised wavefunction of the quantum harmonic oscillator. Since $x$ is the variable representing position in the one-dimensional space that the oscillator is moving in, the normalisation of the wavefunction refers to an integral with respect to $x$ , i.e.

$N^2\int_{-\infty}^{\infty}e^{-\frac{m\omega}{\hbar}x^2}\biggr\[H_n\biggr$\sqrt{\frac{m\omega}{\hbar}}x\biggr$\biggr\]^2dx=1\;\;\;\;\;\;\;\;45$

Using ${dx=\biggr$\frac{\hbar}{m\omega}\biggr$^{\frac{1}{2}}}dy$ , eq45 is equivalent to

$N^2\biggr$\frac{\hbar}{m\omega}\biggr$^{\frac{1}{2}}\int_{-\infty}^{\infty}e^{-y^2}\[H_n(y)\]^2dy=1$

Therefore,

$N^2\biggr$\frac{\hbar}{m\omega}\biggr$^{\frac{1}{2}}2^nn!\sqrt{\pi}=1\;\;or\;\;N=\sqrt[4]{\frac{m\omega}{\pi\hbar}}\frac{1}{\sqrt{2^nn!}}$

The normalised wavefunction is then

$\psi_n(y)=\sqrt[4]{\frac{m\omega}{\pi\hbar}}\frac{1}{\sqrt{2^nn!}}H_n(y)e^{-\frac{y^2}{2}}\;\;\;n=0,1,2,\cdots\;\;\;\;\;\;\;\;45$

$\psi_n(x)=\sqrt[4]{\frac{m\omega}{\pi\hbar}}\frac{1}{\sqrt{2^nn!}}H_n\biggr$\sqrt{\frac{m\omega}{\hbar}}x\biggr$e^{-\frac{m\omega}{2\hbar}x^2}\;\;\;n=0,1,2,\cdots\;\;\;\;\;\;\;\;46$

Question

Show that $\psi_n(x)$ is either an even function or an odd function.

Answer

An even function is defined as $f(x)=f(-x)$ , while and odd function is when $-f(x)=f(-x)$ . So, $e^{-\frac{m\omega}{2\hbar}x^2}$ is an even function. On the other hand, $H_n\biggr$\sqrt{\frac{m\omega}{\hbar}}x\biggr$$ is an even function if $n$ is even and an odd function if $n$ is odd. Since the product of two even functions is an even function, while the product of an odd function and an even function is an odd function, $\psi_n(x)$ is an even function when $n$ is even and an odd function if $n$ is odd.

Next article: Vibration of diatomic molecules

Previous article: Generating function and Rodrigues’ formula for the Hermite polynomials

Content page of vibrational spectroscopy

Content page of advanced chemistry

Main content page

April 30, 2024November 8, 2024

The Born-Oppenheimer approximation

The Born-Oppenheimer approximation assumes that the wavefunctions of atomic nuclei and electrons in a molecule can be treated separately, primarily because the nuclei are much heavier than the electrons.

We begin with the molecular Hamiltonian operator $\hat{H}$ :

$\hat{H}=-\sum_{\alpha}\frac{\hbar^2}{2m_{\alpha}}\nabla_{\alpha}^2-\sum_{i}\frac{\hbar^2}{2m_{e}}\nabla_{i}^2+\hat{V}_{NN}+\hat{V}_{Ne}+\hat{V}_{ee}\;\;\;\;\;\;\;\;51$

where

- $N$ and $e$ refer to nuclei and electrons respectively.
- The first and second terms are the kinetic energy operator of nuclei and electrons respectively.
- $\hat{V}$ is the potential energy of interaction between two particles.

For example, the Hamiltonian for $H_2$ (see diagram above) is:

$\begin{align}\hat{H}=&-\sum_{\alpha}\frac{\hbar^2}{2m_{\alpha}}\nabla_{\alpha}^2-\sum_{\beta}\frac{\hbar^2}{2m_{\beta}}\nabla_{\beta}^2-\frac{\hbar^2}{2m_{e}}\nabla_{1}^2-\frac{\hbar^2}{2m_{e}}\nabla_{2}^2-\frac{e^2}{4\pi\varepsilon_0r_{1\alpha}}\\&-\frac{e^2}{4\pi\varepsilon_0r_{1\beta}}-\frac{e^2}{4\pi\varepsilon_0r_{2\alpha}}-\frac{e^2}{4\pi\varepsilon_0r_{2\beta}}+\frac{e^2}{4\pi\varepsilon_0r_{12}}+\frac{e^2}{4\pi\varepsilon_0R}\end{align}$

To solve the complicated Schrodinger equation $\hat{H}\psi=E\psi$ for a molecule, Max Born and J. Oppenheimer proposed separating the electronic and nuclear motions. The first assumption the two scientists made is that the molecular wavefunction $\psi$ is a product of the wavefunction of nuclear motion $\psi_N(q_{\alpha})$ and the wavefunction of electron motion $\psi_e(q_i,q_{\alpha})$ , i.e.

$\psi=\psi_N\psi_e\;\;\;\;\;\;\;\;52$

Question

Why is $\psi_e$ a function of nuclear coordinates $q_{\alpha}$ and electronic coordinates $q_i$ ?

Answer

The electron density of a molecule changes when the molecule vibrates and therefore the form of the wavefunction of electron motion changes based on the positions of the electrons and the nuclei.

With regard to eq51 and eq52, the Schrodinger equation is

$\biggr$-\sum_{\alpha}\frac{\hbar^2}{2m_{\alpha}}\nabla_{\alpha}^2-\sum_{i}\frac{\hbar^2}{2m_{e}}\nabla_{i}^2+\hat{V}_{NN}+\hat{V}_{Ne}+\hat{V}_{ee}\biggr$\psi_N\psi_e=E\psi_N\psi_e\;\;\;\;\;\;\;\;53$

Using the product rule, we have $\nabla_{\alpha}^2\psi_N\psi_e=\psi_e\nabla_{\alpha}^2\psi_N+2(\nabla_{\alpha}\psi_e)\nabla_{\alpha}\psi_N+\psi_N\nabla^2_{\alpha}\psi_e$ , which when substituted in eq53 gives

$\begin{align}&-\sum_{\alpha}\frac{\hbar^2}{m_{\alpha}}(\nabla_{\alpha}\psi_e)\nabla_{\alpha}\psi_N-\sum_{\alpha}\frac{\hbar^2}{2m_{\alpha}}\psi_N\nabla^2_{\alpha}\psi_e-\sum_{\alpha}\frac{\hbar^2}{2m_{\alpha}}\psi_e\nabla^2_{\alpha}\psi_N\\&-\psi_N\sum_{i}\frac{\hbar^2}{2m_{e}}\nabla_{i}^2\psi_e+(\hat{V}_{NN}+\hat{V}_{Ne}+\hat{V}_{ee})\psi_N\psi_e=E\psi_N\psi_e\;\;\;\;\;\;\;\;54\end{align}$

The second assumption Born and Oppenheimer made was that $\psi_e$ is a slowly varying function of $q_{\alpha}$ . This is because the heavier nuclei are considered relatively stationary as they move much more slowly than electrons. Therefore, the change in the distribution of the electrons with respect to the positions of the nuclei $\nabla_{\alpha}\psi_e$ is very small and eq54 becomes

$-\psi_e\sum_{\alpha}\frac{\hbar^2}{2m_{\alpha}}\nabla^2_{\alpha}\psi_N-\psi_N\sum_{i}\frac{\hbar^2}{2m_{e}}\nabla_{i}^2\psi_e+(\hat{V}_{NN}+\hat{V}_{Ne}+\hat{V}_{ee})\psi_N\psi_e=E\psi_N\psi_e\;\;\;\;\;\;\;\;55$

Substituting the purely electronic Hamiltonian $\hat{H}_e=-\sum_{i}\frac{\hbar^2}{2m_{e}}\nabla_{i}^2+\hat{V}_{Ne}+\hat{V}_{ee}$ and its corresponding Schrodinger equation $\hat{H}_e\psi_e=E_e\psi_e$ in eq55, we have

$\biggr$-\sum_{\alpha}\frac{\hbar^2}{2m_{\alpha}}\nabla^2_{\alpha}+E_e+\hat{V}_{NN}\biggr$\psi_N=E\psi_N\;\;\;\;\;\;\;\;56$

Question

How do we interpret the term $E_e+\hat{V}_{NN}$ in eq56?

Answer

$E_e$ is the eigenvalue of the purely electronic Schrodinger equation $\hat{H}_e\psi_e=E_e\psi_e$ . Using the second assumption, we regard the nuclei as fixed while the electrons move. $q_{\alpha}$ in $\psi_e$ then becomes a constant for a particular nuclear configuation. In other words, $\psi_e$ is now a function of $q_i$ and at the same time parametrically dependent on $q_{\alpha}$ . For each electronic state , the purely electronic Schrodinger equation is solved repeatedly using different fixed values of $q_{\alpha}$ , giving a set of $\psi_e$ and a corresponding set of $E_e$ . Since $\hat{V}_{NN}=\sum_{\alpha}\sum_{\beta >\alpha}\frac{Z_{\alpha}Z_{\beta}e^2}{4\pi\epsilon_0r_{\alpha\beta}$ , we can then easily calculate a value $U=E_e+\hat{V}_{NN}$ for every $q_{\alpha}$ and obtain a curve (see diagram below). In effect, $U$ is a function of $q_{\alpha}$ , making it the potential energy for nuclear motion:

$\biggr\[-\sum_{\alpha}\frac{\hbar^2}{2m_{\alpha}}\nabla^2_{\alpha}+U(q_{\alpha})\biggr\]\psi_N=E\psi_N\;\;\;\;\;\;\;\;57$

We have therefore transformed eq53, the molecular Schrodinger equation, to eq57, which is an eigenvalue equation of a single variable of nuclear coordinates, with the eigenfunction $\psi_N$ as the wavefunction of nuclear motion, the eigenvalue $E$ representing the total energy of the molecule, and the term $U(q_{\alpha})$ as the potential energy for nuclear motion.

Question

How does the above molecular potential energy curve relate to the quadratic potential energy curve derived using the harmonic oscillator model?

Answer

The two curves coincide at low vibration energy levels, i.e., around the equilibrium internuclear distance (see diagram below). Nuclear motion includes translational, rotational and vibrational motions. As translational and rotational motions are purely kinetic, $U(q_{\alpha})$ in eq57 is associated only with molecular vibrational motions. Consider a diatomic molecule with one of the nuclei at the origin and the other on the $x$ -axis. The potential energy term becomes $U(R)$ , where $R$ is the internuclear distance. Using a Taylor series, we can approximate the potential function by expanding it around its minimum at the equilibrium internuclear distance $R=R_e$ :

$U(R)=U(R_e)+U'(R_e)(R-R_e)+\frac{1}{2}U''(R_e)(R-R_e)^2+\cdots$

At $R_e$ , the slope of $U(R)$ is zero, i.e. $U'(R_e)=0$ . For small vibrations, where $R$ is close to $R_e$ , we ignore all the higher power terms, giving $U(R)=U(R_e)+\frac{1}{2}U''(R_e)(R-R_e)^2$ . The physical observables of a system, such as the transition frequencies $\Delta\nu$ of vibrational modes of a molecule, are not affected by the $U(R_e)$ term because it cancels out when we take the difference of two energy levels. Therefore, we can set $U(R_e)$ to zero, resulting in $U(R)=\frac{1}{2}U''(R_e)(R-R_e)^2$ , which is equivalent to the potential energy term of the harmonic oscillator (see eq49), where $U''(R_e)=k$ and $r=R-R_e$ .

In conclusion, the Born-Oppenheimer approximation offers a different approach, compared to the quantum harmonic oscillator, in deriving the Schrodinger equation for a diatomic molecule (see the Hamiltonian of eq57 vs eq48) and the potential energy term for the vibrational motion of a diatomic molecule.

Next article: Nuclear motion in diatomic molecules

Previous article: Vibration of diatomic molecules

Content page of vibrational spectroscopy

Content page of advanced chemistry

April 30, 2024

Vibration of diatomic molecules

The vibration of a diatomic molecule can be modelled using the quantum harmonic oscillator.

The Schrodinger equation can be written as:

$\frac{p_1^2}{2m_1}\psi(x_1,x_2)+\frac{p_2^2}{2m_2}\psi(x_1,x_2)+\frac{1}{2}k(x_2-x_1-x_0)^2\psi(x_1,x_2)=E\psi(x_1,x_2)\;\;\;\;\;\;\;\;47$

To solve the equation, we need to make a change of coordinates to achieve a separation of variables. The new coordinate system consists of the variables $r=x_2-x_1-x_0$ and $s=\frac{m_1x_1+m_2x_2}{m_1+m_2}$ , where $s$ is the centre of mass of the molecule.

Question

What do the terms in eq47 mean? What is a separation of variables?

Answer

Eq47 is an eigenvalue equation, where the Hamiltonian $\hat{H}=\frac{\hat{p}_1^2}{2m_1}+\frac{\hat{p}_2^2}{2m_2}+\frac{1}{2}k(x_2-x_1-x_0)^2$ has two kinetic energy terms, one for each atom, and a potential energy term that describes the interaction between the atoms. According to the classical harmonic oscillator, the potential energy of the system is equal to $\frac{1}{2}kx^2$ , where $x$ is the displacement of the spring from its equilibrium position. For the diatomic molecule in the diagram above, $x=r=x_2-x_1-x_0$ .

To solve eq47, we need it to be of the form $\[\hat{H}_1(x_1)+\hat{H}_2(x_2)\]\psi(x_1,x_2)=(E_1+E_2)\psi(x_1,x_2)$ , where the Hamiltonian is separated into two, each being a function of only one variable. Unless it is transformed, the potential energy term in eq47 does not allow the separation to occur.

The derivatives of $r$ and $s$ with respect to time are $\dot r=\dot x_2-\dot x_1$ and $\dot s=\frac{m_1\dot x_1+m_2\dot x_2}{m_1+m_2}$ respectively. Substituting $\dot x_1=\dot x_2-\dot r$ and $\dot x_x=\dot x_1+\dot r$ in $\dot s$ and rearranging, we have $\dot x_2=\dot s+\frac{m_1\dot r}{m_1+m_2}$ and $\dot x_1=\dot s+\frac{m_2\dot r}{m_1+m_2}$ respectively. Using $p_i=m_i\dot x_i$ , the Hamiltonian of eq47 becomes

$\hat{H}=\frac{1}{2}(M\dot s^2+\mu\dot r^2)+\frac{1}{2}kr^2=\frac{p_s^2}{2M}+\frac{p_r^2}{2\mu}+\frac{1}{2}kr^2\;\;\;\;\;\;\;\;48$

where $p_s=M\dot s$ , $M=m_1+m_2$ , $p_r=\mu\dot r$ and $\mu=\frac{m_1m_2}{m_1+m_2}$ .

Since $M$ and $\dot s$ are the total mass of the molecule and the centre of mass coordinate of the molecule respectively, $\frac{p_s^2}{2M}$ corresponds to the kinetic energy of the translational motion of the molecule. The other kinetic energy energy term $\frac{p_r^2}{2\mu}$ must correspond to the internal motion (rotational and vibrational motion) of the molecule. However, there are only two degrees of freedom for a diatomic molecule in a one-dimensional space – a translational degree of freedom and a vibrational degree of freedom. Discarding the kinetic energy term for the translational motion, the Schrodinger equation describing the vibrational motion of the molecule is

$\biggr$\frac{p_r^2}{2\mu}+\frac{1}{2}kr^2\biggr$\psi(r)=E\psi(r)\;\;\;\;\;\;\;\;49$

Eq49 is equivalent to eq4 except for the substitution of $\mu$ , which is called the reduced mass. Its solutions are given by eq45 (with $m$ replaced with $\mu$ ). With reference to eq21, the vibrational energies of the diatomic molecule are:

$E_n\biggr$n+\frac{1}{2}\biggr$\hbar\omega\;\;\;\;n=0,1,2,\cdots\;\;\;\;\;\;\;\;50$

where ${\omega=\sqrt{\frac{k}{\mu}}}$ .

The potential energy curve of $\frac{1}{2}kr^2$ in eq49 for a typical diatomic molecule is represented by the dashed curve in the diagram below. The minimum potential energy of the molecule occurs when the atoms are at their equilibrium positions. In the next article, we will show that the actual potential energy curve has the shape of the solid curve, which correctly shows that the potential energy of the molecule increases rapidly as the bond length shortens and approaches a constant value representing the dissociation energy of the molecule as $r$ increases. Despite that, the two curves coincide at low vibrational energy levels,. This implies that the harmonic oscillator model of a diatomic molecule only serves as a good approximation at low vibrational energy levels.

Next article: The Born-oppenheimer approximation

Previous article: Normalisation constant for the wavefunctions of the quantum harmonic oscillator

Content page of vibrational spectroscopy

Content page of advanced chemistry

April 30, 2024April 10, 2025

Nuclear motion in diatomic molecules

Nuclear motion in diatomic molecules involves the rotation, vibration and overall translation of nuclei, significantly impacting molecular behaviour and spectroscopic transitions.

From eq57, the Schrodinger equation for nuclear motion in a diatomic-molecule is

$\biggr\[-\frac{\hbar^2}{2m_{\alpha}}\nabla^2_{\alpha}-\frac{\hbar^2}{2m_{\beta}}\nabla^2_{\beta}+U(R)\biggr\]\psi_N=E\psi_N\;\;\;\;\;\;\;\;60$

where $\alpha$ and $\beta$ refer to the nuclei and $R$ is the inter-nuclear distance.

The eigenfunctions and eigenvalues of eq60, which is a two-particle problem, are $E=E_{tr}+E_{int}$ and $\psi_N=\psi_{tr}\psi_{int}$ respectively, where the subscripts tr and int refer to translational motion and internal motion respectively. Since the operator $\frac{\hat{p}^2_{\mu}}{2\mu}$ in eq49 is equivalent to $-\frac{\hbar^2}{2\mu}\nabla^2$ (see this article for derivation), the Schrodinger equation for internal motion is

$\biggr\[-\frac{\hbar^2}{2\mu}\nabla^2+U(R)\biggr\]\psi_{int}=E_{int}\psi_{int}\;\;\;\;\;\;\;\;61$

The potential energy $U$ is spherically symmetric, as it depends on $R$ only. This implies that eq61 is a central force problem, with the following Schrodinger equation and eigenfunction:

$\biggr\[-\frac{\hbar^2}{2\mu}\biggr$\frac{2}{R}\frac{\partial}{\partial R}+\frac{\partial^2}{\partial R^2}\biggr$+\frac{\hat{L}^2}{2\mu R^2}+U(R)\biggr\]\psi_{int}=E_{int}\psi_{int}\;\;\;\;\;\;\;\;62$

$\psi_{int}=P(R)Y^M_J(\theta,\phi)\;\;\;\;\;\;\;\;63$

where $\hat{L}^2$ is the operator for the square of the magnitude of the orbital angular momentum of the molecule, $P(R)$ is a function of $R$ and $Y^M_J(\theta,\phi)$ is a spherical harmonic function with quantum numbers $J$ and $M$ .

The kinetic energy of internal motion is the sum of rotational kinetic energy and vibrational kinetic energy. It follows that the Hamiltonian of eq62 consists, in part, of two kinetic energy operators. If we assume that the diatomic molecule is a rigid rotor ( $R$ is a constant), whose rotational motion is independent of its vibrational motion, then there is no potential energy contribution in the Hamiltonian for the rotation of a rigid rotor because the bond length remains fixed at its equilibrium distance. In other words, rotation is characterised only by $\theta$ and $\phi$ , and the Hamiltonian for vibrational motion consists of both a kinetic energy operator and a potential energy function, as the vibrational motion of the molecule is described by a change in $R$ .

Since ${\biggr\[\frac{\hat{L}^2}{2\mu R^2}Y^M_J(\theta,\phi)\biggr\]P(R)=E_{rot}P(R)Y^M_J(\theta,\phi)}$ , where ${E_{rot}=\frac{J(J+1)\hbar^2}{2\mu R^2}}$ (see this article for derivation)

$\biggr\[-\frac{\hbar^2}{2\mu}\biggr$\frac{2}{R}\frac{\partial}{\partial R}+\frac{\partial^2}{\partial R^2}\biggr$+U(R)\biggr\]P(R)Y^M_J(\theta,\phi)=(E_{int}-E_{rot})P(R)Y^M_J(\theta,\phi)$

$Y^M_J(\theta,\phi)$ becomes a non-zero constant because the Hamiltonian in the above equation is a function of $R$ . So,

$\biggr\[-\frac{\hbar^2}{2\mu}\biggr$\frac{2}{R}\frac{\partial}{\partial R}+\frac{\partial^2}{\partial R^2}\biggr$+U(R)\biggr\]P(R)=(E_{int}-E_{rot})P(R)\;\;\;\;\;\;\;\;64$

To simplify eq64, we substitute ${P(R)=\frac{F(R)}{R}}$ , ${P'(R)=\frac{F'(R)}{R}-\frac{F(R)}{R^2}}$ and ${P''(R)=\frac{F''(R)}{R}-2\frac{F'(R)}{R^2}}+2\frac{F(R)}{R^3}$ in eq64 to give

$\biggr\[-\frac{\hbar^2}{2\mu}\frac{d^2}{dR^2}+U(R)\biggr\]F(R)=(E_{int}-E_{rot})F(R)\;\;\;\;\;\;\;\;65$

We can approximate $U(R)$ by expanding it in a Taylor series around the equilibrium internuclear distance $R=R_e$ :

$U(R)=U(R_e)+U'(R_e)(R-R_e)+\frac{1}{2}U''(R_e)(R-R_e)^2+\cdots\;\;\;\;\;\;\;\;66$

At $R_e$ , the slope of $U(R)$ is zero, i.e. $U'(R_e)=0$ . For small vibrations, where $R$ is close to $R_e$ , we ignore all the higher power terms, giving ${U(R)=U(R_e)+\frac{1}{2}U''(R_e)(R-R_e)^2}$ , which when substituted in eq65 gives

$\biggr\[-\frac{\hbar^2}{2\mu}\frac{d^2}{dR^2}+U(R_e)+\frac{1}{2}U''(R_e)(R-R_e)^2\biggr\]F(R)=(E_{int}-E_{rot})F(R)$

The molecule is in its most stable state when $R=R_e$ for a particular electronic configuration, with $U(R_e)$ being a minimum point on the potential energy curve. This suggests that the electrons are in their lowest energy state for that configuration, and $U(R_e)$ is called the equilibrium electronic energy $E_{elec}$ of the molecule.

Question

Elaborate on the concept of $E_{elec}$ .

Answer

$E_{elec}$ is the energy of the electrons that are interacting with the nuclei at their equilibrium positions. It is neither the purely electronic energy $E_{e}$ nor $E_{e}+\hat{V}_{NN}$ in eq56. It is also not a form of vibrational energy, which is characterised by a displacement of the nuclei from their equilibrium positions. Since $E_{elec}$ is a constant for a particular electronic state, it is subtracted out when we calculate vibrational transition energies for a given electronic state. This is why it is sometimes set to zero.

Letting $x=R-R_e$ and $k=U''(R_e)$ and therefore $F(R)=F(x+R_e)=S(x)$ , we have

$\biggr\[-\frac{\hbar^2}{2\mu}\frac{d^2}{dx^2}+\frac{1}{2}kx^2\biggr\]S(x)=E_{vib}S(x)\;\;\;\;\;\;\;\;67$

where $E_{vib}=E_{int}-E_{rot}-E_{elec}$ .

Question

Show that ${\frac{d^2}{dR^2}F(R)=\frac{d^2}{dx^2}S(x)}$ .

Answer

According to the chain rule, ${\frac{dS(x)}{dx}=\frac{dF(R)}{dx}=\frac{dF(R)}{dR}\frac{dR}{dx}}$ . Since $x=R-R_e$ , we have ${\frac{dR}{dx}=1}$ and ${\frac{dS(x)}{dx}=\frac{dF(R)}{dR}}$ . So, ${\frac{d^2S(x)}{dx^2}=\frac{d}{dx}F'(R)=\frac{dF'(R)}{dR}\frac{dR}{dx}=\frac{d^2F(R)}{dR^2}}$ .

Eq67 is the Schrodinger equation for a one-dimensional harmonic oscillator. Assuming that $S(x)$ behaves the same as a wave function of the harmonic oscillator, it is given by eq46. $E_{vib}$ is then represented by eq50. Therefore, $E$ and $\psi_N$ in eq60 are

$E=E_{tr}+E_{rot}+E_{vib}+E_{elec}\;\;\;\;\;\;\;\;68$

$\psi_N=\psi_{tr}\psi_{rot}\psi_{vib}\;\;\;\;\;\;\;\;69$

where $\psi_{vib}=\frac{S(x)}{R}$ .

Next article: Anharmonicity

Previous article: The Born-Oppenheimer approximation

Content page of vibrational spectroscopy

Content page of advanced chemistry

April 30, 2024October 2, 2025

Vibration of polyatomic molecules

The vibration of polyatomic molecules involves stretching, bending and torsional motions.

The classical kinetic energy $T$ of a vibrating $N$ -nuclei molecule is

$T=\frac{1}{2}\sum_{\alpha=1}^Nm_{\alpha}\biggr\[\biggr$\frac{dx_{\alpha}}{dt}\biggr$^2+\biggr$\frac{dy_{\alpha}}{dt}\biggr$^2+\biggr$\frac{dz_{\alpha}}{dt}\biggr$^2\biggr\]\;\;\;\;\;\;\;\;85$

where $x_{\alpha}$ , $y_{\alpha}$ and $z_{\alpha}$ are the displacements of the $\alpha$ -th nucleus from its equilibrium position of $(x_{\alpha,e},y_{\alpha,e},z_{\alpha,e})$ (see diagram below).

Eq85 can be simplified by changing the Cartesian displacement coordinates $x_{\alpha},y_{\alpha},z_{\alpha}$ to mass-weighted displacement coordinates $q_1,\cdots,q_{3N}$ , where

$\begin{matrix} q_1=m_1^{1/2}x_1&q_2=m_1^{1/2}y_1&q_3=m_1^{1/2}z_1\\q_4=m_2^{1/2}x_2&\cdots &q_{3N}=m_N^{1/2}z_N\end{matrix}$

which gives

$T=\frac{1}{2}\sum_{i=1}^{3N}\dot q^2_i=\frac{1}{2}(\dot q_1\cdots \dot q_{3N})\begin{pmatrix}\dot q_1\\\vdots\\\dot q_{3N}\end{pmatrix}=\frac{1}{2}\boldsymbol{\mathit{\dot q}}^T\boldsymbol{\mathit{\dot q}}\;\;\;\;\;\;\;\;86$

where $\dot q_{i}=\frac{dq_i}{dt}$ and $\boldsymbol{\mathit{\dot q}}$ is a column vector with elements $\dot q_1,\cdots,\dot q_{3N}$ .

The potential energy $U$ of a vibrating $N$ -nuclei molecule is also a function of $x_{\alpha},y_{\alpha},z_{\alpha}$ and hence a function of $q_1,\cdots,q_{3N}$ . We can approximate the multivariable function $U$ by expanding it in a Maclaurin series around the equilibrium positions:

$U(q_1,\cdots,q_{3N})=U(q_{1,e=0},\cdots,q_{3N,e=0})+\sum_{i=1}^{3N}\frac{\partial U(q_{1,e=0},\cdots,q_{3N,e=0})}{\partial q_i}q_i+\frac{1}{2!}\sum_{i,j=1}^{3N}\frac{\partial^2 U(q_{1,e=0},\cdots,q_{3N,e=0})}{\partial q_i\partial q_j}q_iq_j+\cdots$

which can be abbreviated to

$U=U_e+\sum_{i=1}^{3N}\biggr$\frac{\partial U}{\partial q_i}\biggr$_eq_i+\frac{1}{2!}\sum_{i,j=1}^{3N}\biggr$\frac{\partial^2 U}{\partial q_i\partial q_j}\biggr$_eq_iq_j+\cdots$

When the nuclei are at their equilibrium positions, the slope of $U$ is zero, i.e. ${\biggr$\frac{\partial U}{\partial q_i}\biggr$_e=0}$ . For small vibrations, where $q_i$ is close to $q_{i,e}$ , we ignore all the higher power terms, giving

$U=U_e+\frac{1}{2}\sum_{i,j=1}^{3N}\biggr$\frac{\partial^2 U}{\partial q_i\partial q_j}\biggr$_eq_iq_j\;\;\;\;\;\;\;\;87$

or in matrix notation

$U=U_e+\frac{1}{2}\boldsymbol{\mathit{q}}^T\boldsymbol{\mathit{Uq}}\;\;\;\;\;\;\;\;88$

where $\boldsymbol{\mathit{q}}$ is the column vector with elements $q_1,\cdots,q_{3N}$ and $\boldsymbol{\mathit{U}}$ is a $3N\times 3N$ matrix, called the Hessian matrix, with elements $u_{ij}=\small{\biggr$\frac{\partial^2 U}{\partial q_i\partial q_j}\biggr$_e}$ .

Question

Show that $\boldsymbol{\mathit{q}}^T\boldsymbol{\mathit{Uq}}=\sum_{i,j=1}^{3N}\biggr$\frac{\partial^2 U}{\partial q_i\partial q_j}\biggr$_eq_iq_j$ .

Answer

$\begin{pmatrix}q_1&q_2&\cdots&q_{3N}\end{pmatrix}\begin{pmatrix}u_{11}&u_{12}&\cdots&u_{1j}\\u_{21}&u_{22}&\cdots&u_{2j}\\\vdots&\vdots&\ddots&\vdots\\u_{i1}&u_{i2}&\cdots&u_{ij}\end{pmatrix}\begin{pmatrix}q_1\\q_2\\\vdots\\q_{3N}\end{pmatrix}=\\\begin{pmatrix}q_1&q_2&\cdots&q_{3N}\end{pmatrix}\begin{pmatrix}\sum_{j=1}^{3N}u_{1j}q_j\\\sum_{j=1}^{3N}u_{2j}q_j\\\vdots\\\sum_{j=1}^{3N}u_{ij}q_j\end{pmatrix}=\sum_{i,j=1}^{3N}u_{ij}q_iq_j$

We want to derive an eigenvalue equation that resembles the harmonic oscillator equation so that we can easily solve for the vibrational energies of a polyatomic molecule. However, the Hamiltonian formed by combining eq86 and eq87 doesn’t resemble the Hamiltonian of a harmonic oscillator. Therefore, we need to further simply eq86 and eq87 with a change of variables.

With regard to eq88, the potential energy of a system is a real-valued physical property. This implies that $\boldsymbol{\mathit{U}}$ must consists of real-valued entries $u_{ij}$ . Since $u_{ij}=\small{\biggr$\frac{\partial^2 U}{\partial q_i\partial q_j}\biggr$_e=\biggr$\frac{\partial^2 U}{\partial q_j\partial q_i}\biggr$_e=u_{ji}}$ , the matrix $\boldsymbol{\mathit{U}}$ is both real and symmetric.

Question

i) Prove that eigenvectors of a real and symmetric matrix $O$ with distinct eigenvalues are orthogonal.

ii) What if some of the eigenvalues are degenerate?

Answer

i) Let $Ox=ax$ and $Oy=by$ , where $a\neq b$ . Multiply the first eigenvalue equation by $y^T$ , we have

$y^Tax=(y^TO)x=(O^Ty)^Tx=by^Tx$

where the 2^nd equality uses the identity mentioned in this article, while the 3^rd equality employs the property of a symmetrc matrix.

Since $a\neq b$ , we have $y^Tx=0$ , i.e. $y$ is orthogonal to $x$ .

ii) When some eigenvalues are degenerate, the corresponding eigenvectors belong to the same eigenspace. Since we can always select a set of linearly independent eigenvectors that span the eigenspace, these eigenvectors can be orthogonalised using the Gram-Schmidt process without changing the eigenvalue. So, for a real symmetric matrix, we can always find a basis of eigenvectors that are orthogonal to each other, regardless of whether the eigenvalues are distinct or degenerate.

Let the orthogonal eigenvectors $\boldsymbol{\mathit{l}}_m$ of $\boldsymbol{\mathit{U}}$ be the columns of the matrix $\boldsymbol{\mathit{L}}$ , where $\boldsymbol{\mathit{L}}^T\boldsymbol{\mathit{L=I}}$ . As shown in this article, we can express the eigenvalue problem $\boldsymbol{\mathit{Ul}}_m=\lambda_j\boldsymbol{\mathit{l}}_m$ as

$\boldsymbol{\mathit{UL=L\Lambda}}\;\;\;or\;\;\;\boldsymbol{\mathit{L}}^T\boldsymbol{\mathit{UL=\Lambda}}$

where $\boldsymbol{\mathit{\Lambda}}$ is the diagonal eigenvalue matrix with diagonal entries $\lambda_m$ .

Computing the secular equation $det(u_{ij}-\lambda_m\delta_{ij})$ allows us to find $\lambda_m$ , and solving the simultaneous equations of $((\boldsymbol{\mathit{U}}-\lambda_m \boldsymbol{\mathit{I}})\boldsymbol{\mathit{l}}_m=0$ , which is equivalent to $\sum_{j=1}^{3N}(u_{ij}-\lambda_m\delta_{ij})l_{jm}=0$ (where $i=1,\cdots,3N$ , enables us to find $\boldsymbol{\mathit{l}}_m$ .

Consider the following linear combinations of mass-weighted displacement coordinates

$Q_k=\sum_{j=1}^{3N}l_{jk}q_{j}=\boldsymbol{\mathit{L}}^T\boldsymbol{\mathit{q=Q}}\;\;\;\;\;\;k=1,\cdots,3N\;\;\;\;\;\;\;\;89$

where $\boldsymbol{\mathit{Q}}$ is the column vector with elements $Q_1,\cdots,Q_{3N}$ .

Since $l_{jk}$ are the components of the $k$ -th orthogonal eigenvector of $\boldsymbol{\mathit{U}}$ and each basis $q_j$ describes the displacement of a nucleus, $Q_k$ , known as the normal coordinates, is a set of orthogonal vectors that expresses $3N$ independent motions (called degrees of freedom) of the molecule as a whole. In other words, the normal coordinates measure the motion of atoms in a molecule, where all the atoms move concertedly with respect to their equilibrium positions.

Substituting $\boldsymbol{\mathit{q=LQ}}$ in eq88

$U=U_e+\frac{1}{2}(\boldsymbol{\mathit{LQ}})^T\boldsymbol{\mathit{ULQ}}=U_e+\frac{1}{2}\boldsymbol{\mathit{Q}}^T\boldsymbol{\mathit{L}}^T\boldsymbol{\mathit{ULQ}}=U_e+\frac{1}{2}\boldsymbol{\mathit{Q}}^T\boldsymbol{\mathit{\Lambda Q}}=U_e+\frac{1}{2}\sum_{j=1}^{3N}\lambda_iQ_j^2\;\;\;\;\;\;\;\;90$

Since $\boldsymbol{\mathit{q=LQ}}=\sum_{j=1}^{3N}l_{kj}Q_j=q_k$ , we have $\frac{dq_k}{dt}=\sum_{j=1}^{3N}l_{kj}\frac{dQ_j}{dt}$ or $\boldsymbol{\mathit{\dot q=L\dot Q}}$ and eq86 becomes

$T=\frac{1}{2}\boldsymbol{\mathit{\dot q}}^T\boldsymbol{\mathit{\dot q}}=\frac{1}{2}(\boldsymbol{\mathit{L\dot Q}})^T\boldsymbol{\mathit{L\dot Q}}=\frac{1}{2}\boldsymbol{\mathit{\dot Q}}^T\boldsymbol{\mathit{L}}^T\boldsymbol{\mathit{L\dot Q}}=\frac{1}{2}\boldsymbol{\mathit{\dot Q}}^T\boldsymbol{\mathit{\dot Q}}=\frac{1}{2}\sum_{j=1}^{3N}\dot Q_j^2\;\;\;\;\;\;\;\;91$

Combining eq90 and eq91, the classical total energy of the system is $E=\frac{1}{2}\sum_{j=1}^{3N}\dot Q_j^2+\frac{1}{2}\sum_{j=1}^{3N}\lambda_j Q_j^2$ , where we have set $U_e=0$ . This is because we are usually interested in computing transition energies and the eigenvalue corresponding to $U_e$ is subtracted out when we calculate vibrational transition energies for a given electronic state. The Hamiltonian is derived by replacing the classical observables by quantum-mechanical operators:

$\hat{H}=-\frac{\hbar^2}{2}\sum_{j=1}^{3N}\frac{\partial^2}{\partial Q^2_j}+\frac{1}{2}\sum_{j=1}^{3N}\lambda_j Q_j^2\;\;\;\;\;\;\;\;92$

Question

Show that the classical kinetic energy term $\frac{1}{2}\sum_{k=1}^{3N}\dot Q_k^2$ can be replaced by the quantum-mechanical operator $-\frac{\hbar^2}{2}\sum_{k=1}^{3N}\frac{\partial^2}{\partial Q^2_k}$ .

Answer

From eq89, $\dot Q_k=\sum_{j=1}^{3N}l_{kj}\dot q_j$ . Substituting the momentum operator $\hat{p}_1=m_1\dot x_1=\frac{\hbar}{i}\frac{\partial}{\partial x_1}$ in $\dot q_1=m_1^{1/2}\dot x_1$ , we have $\dot q_1=m_1^{-1/2}\frac{\hbar}{i}\frac{\partial}{\partial x_1}$ . Using the single-variable chain rule $\frac{\partial}{\partial x_1}=\frac{\partial}{\partial q_1}\frac{\partial q_1}{\partial x_1}$ , we have $\frac{\partial}{\partial x_1}=m_1^{1/2}\frac{\partial}{\partial q_1}$ , which when substituted in $\dot q_1$ gives $\dot q_1=\frac{\hbar}{i}\frac{\partial }{\partial q_1}$ . We can find similar expressions for $\dot q_2,\cdots,\dot q_{3N}$ . So, $\dot Q_k=\frac{\hbar}{i}\sum_{j=1}^{3N}l_{jk}\frac{\partial}{\partial q_j}$ . Since $U$ is a function of $Q_k$ , and $Q_k$ is a function of $q_j$ , the multi-variable chain rule gives $\frac{\partial}{\partial Q_k}=\sum_{j=1}^{3N}\frac{\partial}{\partial q_j}\frac{\partial q_j}{\partial Q_k}$ . Relabelling $q_k=\sum_{j=1}^{3N}l_{kj}Q_j$ as $q_j=\sum_{k=1}^{3N}l_{jk}Q_k$ , we have $\frac{\partial q_j}{\partial Q_k}=l_{jk}$ and therefore, $\frac{\partial}{\partial Q_k}=\sum_{j=1}^{3N}l_{jk}\frac{\partial}{\partial q_j}$ , which when substituted in $\dot Q_k$ gives $\dot Q_k=\frac{\hbar}{i}\frac{\partial}{\partial Q_k}$ . Substituting this last expression of $\dot Q_k$ in $\frac{1}{2}\sum_{k=1}^{3N}\dot Q_k^2$ gives $-\frac{\hbar^2}{2}\sum_{k=1}^{3N}\frac{\partial^2}{\partial Q_k^2}$ .

The kinetic energy term in eq92 involves a change in $Q_j$ , which is a function of $q_k$ and hence a function of the displacement coordinates $x_{\alpha}$ , $y_{\alpha}$ and $z_{\alpha}$ , which are relative to the equilibrium positions of the nuclei. In other words, $Q_j$ is a function of the displacements of all the nuclei from their respective equilibrium positions. For the molecule as a whole, all the equilibrium positions of the nuclei can be taken as the center of mass of the molecule. Since translational and rotational motions of the molecule involve no displacement relative to the molecule’s centre of mass, $\frac{\partial^2}{\partial Q_j^2}=0$ for such motions. Furthermore, translational and rotational motions are purely kinetic with no potential energy component. Therefore, we can rewrite eq92 as:

$\hat{H}=-\frac{\hbar^2}{2}\sum_{j=1}^m\frac{\partial^2}{\partial Q_j^2}+\frac{1}{2}\sum_{j=1}^m\lambda_jQ_j^2\;\;\;\;\;\;\;\;92$

where $m=3N-6$ and $m=3N-5$ for a non-linear molecule and a linear molecule respectively.

Each $m$ is called a normal mode and corresponds to a unique vibrational motion of the molecule. The Schrodinger equation is

$-\frac{\hbar^2}{2}\sum_{j=1}^m\frac{\partial^2}{\partial Q_j^2}\psi_{vib}+\frac{1}{2}\sum_{j=1}^m\lambda_jQ_j^2\psi_{vib}=E\psi_{vib}\;\;\;\;\;\;\;\;93$

where the separation of variables technique allows us to approximate the total vibrational wavefunction as $\psi_{vib}=\prod_{j=1}^m\psi_j(Q_j)$ and hence, $E=\sum_{j=1}^mE_j$ .

It follows that each $\psi_j(Q_j)$ describes a normal mode of the molecule. Comparing with eq4a, eq93 is the Schrodinger equation for a harmonic oscillator of unit mass and force constant $\lambda_j$ . This implies that (see eq46)

$\psi_j(Q_j)=N_jH_j\biggr$\sqrt{\frac{\omega_j}{\hbar}}Q_j\biggr$e^{-\frac{\omega_jQ^2_j}{2\hbar}}\;\;\;\;\;\;\;\;94$

where $N_j=\sqrt[4]{\frac{\omega_j}{\pi\hbar}}\frac{1}{\sqrt{2^nn!}}$ is the normalisation constant for the Hermite polynomials $H_j$ .

Therefore, eq93 can be solved by computing $m$ individual Schrodinger equations. With reference to eq21, each solution gives

$E_j=\biggr$n_j+\frac{1}{2}\biggr$\hbar\omega_j\;\;\;\;\;\;n_j=0,1,2,\cdots\;\;\;\;\;\;\;\;95$

where $n_j$ is the vibrational quantum number of the $j$ -th normal mode of the molecule.

The total vibrational energy is

$E=\sum_{j=1}^m\biggr$n_j+\frac{1}{2}\biggr$\hbar\omega_j\;\;\;\;\;\;\;\;96$

The vibrational state of a polyatomic molecule $\vert n_1,n_2,\cdots,n_m\rangle$ is characterised by $m$ quantum numbers. For example the vibrational ground state of $H_2O$ , where $m=3N-6=3$ , is $\vert 0,0,0\rangle$ . The ground state vibrational wavefunction of any polyatomic molecule is

$\vert 0,0,\cdots,0\rangle=\psi_0=N\prod_{j=1}^mH_j\biggr$\sqrt{\frac{\omega_j}{\hbar}}Q_j\biggr$e^{-\frac{\omega_jQ_j^2}{2\hbar}}=N\prod_{j=1}^me^{-\frac{\omega_jQ_j^2}{2\hbar}}\;\;\;\;\;\;\;\;97$

Eq97 is a Gaussian function that is symmetrical in $Q_j$ . The implication of this in group theory is that the ground state vibrational wavefunction of any polyatomic molecule is totally symmetric under all symmetry operations of the point group that the molecule belongs to. In fact, group theory provides an easy way to determine the vibrational modes of a molecule and to analyse the possible transitions between different vibrational states.

Next article: Vibrational selection rules for molecules

Previous article: Anharmonicity

Content page of vibrational spectroscopy

Content page of advanced chemistry

April 30, 2024

Recurrence relations of the Hermite polynomials

The recurrence relations of the Hermite polynomials describe how each polynomial in the sequence can be obtained from its predecessors.

Expanding eq18 and assuming $n$ is an even number,

${H_n(y)=(-1)^{\frac{n}{2}}\frac{n!}{(\frac{n}{2})!}\biggr\[1-\frac{2n}{2!}y^2+\frac{2^2n(n-2)}{4!}y^4-\frac{2^3n(n-2)(n-4)}{6!}y^6+\cdots\biggr\]\;\;\;\;\;\;\;\;26}$

${H_{n-1}(y)=(-1)^{\frac{n-2}{2}}\frac{2(n-1)!}{(\frac{n-2}{2})!}\biggr\[y-\frac{2(n-2)}{3!}y^3+\frac{2^2(n-2)(n-4)}{5!}y^5-\frac{2^3(n-2)(n-4)(n-6)}{7!}y^7+\cdots\biggr\]\;\;\;\;\;\;\;\;27}$

Taking the derivative of eq26, we have

${\frac{dH_n(y)}{dy}=(-1)^{\frac{n}{2}-1}\frac{2(n-1)!}{(\frac{n-2}{2})!}\biggr\[2ny-\frac{2^2n(n-2)}{3!}y^3+\frac{2^3n(n-2)(n-4)}{5!}y^5-\cdots\biggr\]\;\;\;\;\;\;\;\;28}$

Substituting eq27 in eq28,

$\frac{dH_n(y)}{dy}=2nH_{n-1}(y)\;\;\;\;\;\;\;\;29$

If we repeat the steps from eq26 through eq29 on the assumption that $n$ is an odd number, we too end up with eq29. Therefore, $n$ in eq29 represents any number. Taking the derivative of eq29 again,

$\frac{d^2H_n(y)}{dy^2}=2n\frac{dH_{n-1}(y)}{dy}=4n(n-1)H_{n-2}(y)\;\;\;\;\;\;\;\;30$
Substitute 29 and 30 in eq23, we have

$H_n(y)-2yH_{n-1}(y)+2(n-1)H_{n-2}(y)=0\;\;\;\;\;\;\;\;31$

Replacing the dummy index $n$ in eq31 with $n+1$ gives

$H_{n+1}(y)-2yH_{n}(y)+2nH_{n-1}(y)=0\;\;\;\;\;\;\;\;32$

Eq29 and eq32 are the recurrence relations of the Hermite polynomials.

Question

Given $H_0(y)=1$ , show that eq32 can be used to generate the Hermite polynomials.

Answer

Substituting $n=0$ in eq32, we have $H_1(y)=2yH_0(y)=2y$ . Substituting $n=1$ in eq32, we have $H_2(y)=4y^2-2$ . Repeating the logic, we can generate the rest of the polynomials.

Question

Using eq32, ${y=\sqrt{\frac{m\omega}{\hbar}}x}$ and eq46, show that ${\langle\psi_k\vert x\vert\psi_n\rangle=\sqrt{\frac{\hbar(n+1)}{2m\omega}}\delta_{k,n+1}+\sqrt{\frac{n\hbar}{2m\omega}}\delta_{k,n-1}}$ .

Answer

Substituting eq32 in eq46, we have

$x\psi_n(x)=\sqrt{\frac{\hbar(n+1)}{2m\omega}}\psi_{n+1}(x)+\sqrt{\frac{n\hbar}{2m\omega}}\psi_{n-1}(x)$

Multiplying the above equation by $\psi^*_k(x)$ and integrating over all space,

$\langle\psi_k\vert x\vert\psi_n\rangle=\sqrt{\frac{\hbar(n+1)}{2m\omega}}\delta_{k,n+1}+\sqrt{\frac{n\hbar}{2m\omega}}\delta_{k,n-1}\;\;\;\;\;\;\;\;32a$

Next article: Taylor series of single-variable and multi-variable functions

Previous article: Orthogonality of the wavefunctions of The quantum harmonic oscillator

Content page of vibrational spectroscopy

Content page of advanced chemistry

Main content page

April 30, 2024August 17, 2025

Anharmonicity

Anharmonicity refers to the deviation of a system’s vibrational motion from the harmonic oscillator model, resulting in changes to the energy levels.

According to the harmonic oscillator model, the potential energy curve of a diatomic molecule is a quadratic function and there are infinite vibrational energy levels that are equally spaced. In reality, the potential energy curve resembles the one derived using the Born-Oppenheimer approximation (see diagram below). Furthermore, the energy levels are finite and the spacing between adjacent levels decreases as $n$ increases.

To understand how anharmonicity results in the converging energy levels, we refer to eq67, which was derived using the Born-Oppenheimer approximation by neglecting some terms in eq66. The inclusion of these terms via the perturbation theory will show the changes in the energy levels.

Consider the following:

Description	Formula	Reference equation
Zeroth-order wavefunction	$\psi^{(0)}_{nJM}=\frac{S_n(x)}{R}Y_J^M(\theta,\phi)$	70
Zeroth-order energy	$E_n^{(0)}=\biggr$n+\frac{1}{2}\biggr$\hbar\omega$	50
1^st-order perturbation Hamiltonian	$\hat{H}^{(1)}=\frac{1}{6}U'''(R_e)x^3+\frac{1}{24}U^{iv}(R_e)x^4$	71
1^st-order perturbation energy	$E_n^{(1)}=\langle\psi_{nJM}^{(0)}\vert\hat{H}^{(1)}\vert\psi_{nJM}^{(0)}\rangle$	72

where $n,J,M$ are quantum numbers and $x=R-R_e$ .

Substituting eq70 in eq72 and noting that $S(x)$ is real,

$E_n^{(1)}=\int_0^{\infty}\frac{S_n(x)^2}{R^2}\biggr\[\frac{1}{6}U'''(R_e)x^3+\frac{1}{24}U^{iv}(R_e)x^4\biggr\]R^2dR\int_0^{2\pi}\int_0^{\pi}[Y_J^M(\theta,\phi)]^2sin\theta d\theta d\phi$

For a set of orthonormal spherical harmonic functions, ${\int_0^{2\pi}\int_0^{\pi}[Y_J^M(\theta,\phi)]^2sin\theta d\theta d\phi=1}$ . Furthermore, ${\frac{dR}{dx}=\frac{d}{dx}(x+R_e)=1\;\;\Rightarrow\;\;dR=dx}$ and so,

$E_n^{(1)}=\int_{-R_e}^{\infty}S_n(x)^2\biggr\[\frac{1}{6}U'''(R_e)x^3+\frac{1}{24}U^{iv}(R_e)x^4\biggr\]dx\;\;\;\;\;\;\;\;73$

If we assume that the probability of a nucleus of a vibrating molecule undergoing a displacement of more than $\vert R_e\vert$ is very small, then $S(x)\approx 0$ for $x < -R_e$ , and we can express the lower limit of the integral in eq73 as $-\infty$ without any major error:

$E_n^{(1)}=\frac{U'''(R_e)}{6}\int_{-\infty}^{\infty}S_n(x)^2x^3dx+\frac{U^{iv}(R_e)}{24}\int_{-\infty}^{\infty}S_n(x)^2x^4dx\;\;\;\;\;\;\;\;74$

As mentioned in this article, $S(x)$ is either an even function or an odd funcion, which makes $S(x)^2$ an even function. Therefore, the integrand $S(x)^2x^3$ is an odd function, whose integral over all space is zero. Eq74 becomes

$E_n^{(1)}=\frac{U^{iv}(R_e)}{24}\langle S_n\vert x^4\vert S_n\rangle\;\;\;\;\;\;\;\;75$

Question

$R$ , $S$ , and $T$ are linear operators with matrix representations $\boldsymbol{\mathit{R}}$ , $\boldsymbol{\mathit{S}}$ and $\boldsymbol{\mathit{T}}$ respectively. The matrix elements of $\boldsymbol{\mathit{R}}$ , $\boldsymbol{\mathit{S}}$ and $\boldsymbol{\mathit{T}}$ are $R_{mn}=\langle f_m\vert R\vert f_n\rangle$ , $S_{mn}=\langle f_m\vert S\vert f_n\rangle$ and $T_{mn}=\langle f_m\vert T\vert f_n\rangle$ respectively, where $\{f_i\}$ is a complete orthonormal basis set. Show that $\boldsymbol{\mathit{R}}=\boldsymbol{\mathit{S}}\boldsymbol{\mathit{T}}$ if $R=ST$ .

Answer

Since $\{f_i\}$ is a complete orthonormal basis set, ${\sum_i\vert f_i\rangle\langle f_i\vert=I}$ , where $I$ is the identity matrix. So,

$\begin{align}R_{mn}&=\langle f_m\vert R\vert f_n\rangle=\langle f_m\vert ST\vert f_n\rangle=\langle f_m\vert S\sum_i\vert f_i\rangle\langle f_i\vert T\vert f_n\rangle\\&=\sum_i\langle f_m\vert S\vert f_i\rangle\langle f_i\vert T\vert f_n\rangle=\sum_iS_{mi}T_{in}\end{align}$

which is the definition of the matrix multiplication of $\boldsymbol{\mathit{R}}=\boldsymbol{\mathit{S}}\boldsymbol{\mathit{T}}$ .

Consider the integral $\langle S_{n'}\vert x^2\vert S_n\rangle$ . As $x^2=xx$ , we have

$\langle S_{n'}\vert x^2\vert S_n\rangle=\sum_i\langle S_{n'}\vert x\vert S_i\rangle\langle S_{i}\vert x\vert S_n\rangle\;\;\;\;\;\;\;\;76$

Replacing $m$ in eq32a with $\mu$ and substituting the resultant equation in eq76

$\begin{align}\langle S_{n'}\vert x^2\vert S_n\rangle=&\frac{\hbar}{2\mu\omega}\biggr\[\sum_i\sqrt{(i+1)(n+1)}\delta_{n',i+1}\delta_{i,n+1}+\sum_i\sqrt{n(i+1)}\delta_{n',i+1}\delta_{i,n-1}\\&+\sum_i\sqrt{i(n+1)}\delta_{n',i-1}\delta_{i,n+1}+\sum_i\sqrt{in}\delta_{n',i-1}\delta_{i,n-1}\biggr\]\;\;\;\;\;\;\;\;77\end{align}$

All the terms in the first summation vanish except for when $i=n'-1$ and $i=n+1$ . Using $i=n+1$ , the 1^st summation becomes $\sqrt{(n+2)(n+1)}\delta_{n',n+2}$ . Applying the same logic to the remaining summations, eq77 becomes

$\langle S_{n'}\vert x^2\vert S_n\rangle=\frac{\hbar}{2\mu\omega}\biggr\[\sqrt{(n+2)(n+1)}\delta_{n',n+2}+(2n+1)\delta_{n',n}+\sqrt{n(n-1)}\delta_{n',n-2}\biggr\]\;\;\;\;\;\;\;\;78$

Let’s now consider the integral $\langle S_{n'}\vert x^4\vert S_n\rangle$ . As $x^4=x^2x^2$ , we have

$\langle S_{n'}\vert x^4\vert S_n\rangle=\sum_i\langle S_{n'}\vert x^2\vert S_i\rangle\langle S_{i}\vert x^2\vert S_n\rangle\;\;\;\;\;\;\;\;79$

Substituting eq78 in eq79 and employing the same logic as before, we have

$\begin{align}\langle S_{n'}\vert x^4\vert S_n\rangle=&\biggr$\frac{\hbar}{2\mu\omega}\biggr$^2\biggr\[\sqrt{(n+4)(n+3)(n+2)(n+1)}\delta_{n',n+4}\\&+(2n+1)\sqrt{(n+2)(n+1)}\delta_{n',n+2}+n(n-1)\delta_{n',n}\\&+(2n+5)\sqrt{(n+2)(n+1)}\delta_{n',n+2}+(2n+1)^2\delta_{n',n}\\&+(2n-3)\sqrt{n(n-1)}\delta_{n',n-2}+(n+2)(n+1)\delta_{n',n}\\&+(2n+1)\sqrt{n(n-1)}\delta_{n',n-2}+\sqrt{n(n-1)(n-2)(n-3)}\delta_{n',n-4}\biggr\]\end{align}$

If $n'=n$ , the above equation becomes

$\langle S_{n}\vert x^4\vert S_n\rangle=\biggr$\frac{\hbar}{2\mu\omega}\biggr$^2(6n^2+6n+3)\;\;\;\;\;\;\;\;80$

Substitute eq80 in eq75,

$E_n^{(1)}=\frac{U^{iv}(R_e)}{8}\biggr$\frac{\hbar}{2\mu\omega}\biggr$^2(2n^2+2n+1)\;\;\;\;\;\;\;\;81$

With reference to eq74, the $x^3$ term disappeared for the computation of the first-order energy correction $E_n^{(1)}$ . So, we need to account for it by calculating the second-order energy correction $E_n^{(2)}$ . From eq270a,

$E_n^{(2)}=\sum_{m\neq n}\frac{\vert\langle S_m^{(0)}\vert\hat{H}^{(1)}\vert S_n^{(0)}\rangle\vert^2}{E_n^{(0)}-E_m^{(0)}}\;\;\;\;\;\;\;\;82$

Substituting eq50 and $H^{(1)}=\frac{1}{6}U'''(R_e)x^3$ in eq82

$E_n^{(2)}=\frac{U'''(R_e)}{6}\sum_{m\neq n}\frac{\vert\langle S_m^{(0)}\vert x^3\vert S_n^{(0)}\rangle\vert^2}{(n-m)\hbar\omega}$

As $x^3=x^2x$ , we have $\langle S_{m}^{(0)}\vert x^3\vert S_n^{(0)}\rangle=\sum_i\langle S_{m}\vert x^2\vert S_i\rangle\langle S_{i}\vert x\vert S_n\rangle$ and

$E_n^{(2)}=\frac{U'''(R_e)}{6}\sum_{m\neq n}\frac{\vert\sum_i\langle S_{m}\vert x^2\vert S_i\rangle\langle S_{i}\vert x\vert S_n\rangle\vert^2}{(n-m)\hbar\omega}\;\;\;\;\;\;\;\;83$

Substituting eq32a (where $m=\mu$ ) and eq78 in eq83, and using the same logic in deriving eq78,

$E_n^{(2)}=\frac{U'''(R_e)}{6}\sum_{m\neq n}\frac{\vert A\vert^2}{(n-m)\hbar\omega}\;\;\;\;\;\;\;\;83$

where

$\begin{align}A=&\frac{\hbar}{2\mu\omega}\sqrt{\frac{(n+1)(n+2)(n+3)\hbar}{2\mu\omega}}\delta_{m,n+3}+3\biggr\[\frac{(n+1)\hbar}{2\mu\omega}\biggr\]^{\frac{3}{2}}\delta_{m,n+1}\\&+3\biggr$\frac{n\hbar}{2\mu\omega}\biggr$^{\frac{3}{2}}\delta_{m,n-1}+\frac{\hbar}{2\mu\omega}\sqrt{\frac{n(n-1)(n-2)\hbar}{2\mu\omega}}\delta_{m,n-3}\end{align}$

Since $\delta_{x,y}\delta_{x,y}=\delta_{x,y}$ and $\delta_{x,y}\delta_{x,z}=0$

$\begin{align}\vert A\vert^2 = & (n+1)(n+2)(n+3)\biggr$\frac{\hbar}{2\mu\omega}\biggr$^3\delta_{m,n+3}+9\biggr\[\frac{(n+1)\hbar}{2\mu\omega}\biggr\]^3\delta_{m,n+1}\\&+9\biggr$\frac{n\hbar}{2\mu\omega}\biggr$^3\delta_{m,n-1}+n(n-1)(n-2)\biggr$\frac{\hbar}{2\mu\omega}\biggr$^3\delta_{m,n-3}\end{align}$

Substituting the above equation in eq83 and expanding the summation gives

$E_n^{(2)}=-\frac{5\hbar^2U'''(R_e)}{8\mu^3\omega^3}\biggr$n^2+n+\frac{11}{30}\biggr$\;\;\;\;\;\;\;\;84$

With reference to eq266, $E_n\approx E_n^{(0)}+\lambda E_n^{(1)}+\lambda^2E_n^{(2)}$ . Since $E_n^{(2)}$ is negative, $\Delta E_n$ becomes smaller as $n$ increases.

Question

Answer

Next article: Generating function and Rodrigues’ formula for the Hermite polynomials

Previous article: Recurrence relations of the Hermite polynomials

Question

Answer

Question

Answer

Question

Answer

Question

Answer

Question

Answer

Question

Answer

Next article: Orthogonality of the wavefunctions of the quantum harmonic oscillator

Previous article: The harmonic oscillator

Question

Answer

Next article: Recurrence relations of the Hermite polynomials

Previous article: The quantum harmonic oscillator

Question

Answer

Question

Answer

Next article: Vibration of diatomic molecules

Previous article: Generating function and Rodrigues’ formula for the Hermite polynomials

Question

Answer

Question

Answer

Question

Answer

Next article: Nuclear motion in diatomic molecules

Previous article: Vibration of diatomic molecules

Question

Answer

Next article: The Born-oppenheimer approximation

Previous article: Normalisation constant for the wavefunctions of the quantum harmonic oscillator

Question

Answer

Question

Answer

Next article: Anharmonicity

Previous article: The Born-Oppenheimer approximation

Question

Answer

Question

Answer

Question

Answer

Next article: Vibrational selection rules for molecules

Previous article: Anharmonicity

Question

Answer

Question

Answer

Next article: Taylor series of single-variable and multi-variable functions

Previous article: Orthogonality of the wavefunctions of The quantum harmonic oscillator

Question

Answer

Next article: Vibration of polyatomic molecules

Previous article: Nuclear motion in diatomic molecules