How is the vibration of multiple atoms treated in vibrational spectroscopy?

Most of the advanced undergraduate textbooks on physical chemistry/ instrumental analysis show the case of a diatomic harmonic oscillator and derive formulae for predicting the frequencies. The concept of reduced mass and force constant is also illustrated for simple atoms. The concept of reduced mass is limited to a two-body problem, although Wikipedia says that it can be generalized for multi-body problems. Reduced mass,

Similarly, the force constant is limited to two bodies; e.g., it is a typical teaching exercise to calculate a force constant from experimental vibrational frequencies.

I was wondering about more complex examples, such as aromatic ring breathing picture from here. The number 1 and 15 is the notation of Wilson for normal modes. For example, how do physical chemists or theoreticians predict the vibrational frequencies where six atoms are executing harmonic motion (i.e., what methods are used to solve such problems?). I wanted to see an elementary exposition of such treatment, perhaps as a solved example in a good reference which does not need computer programs or fancy tools.

Is the generalized reduced mass still useful in polyatomic vibrations, and how is the concept of force constant generalized to more than two atoms? Thanks.

• Well, you're in luck. I happen to have done my doctoral dissertation on IVR (intramolecular vibrational relaxation) for benzene. I will happily send you reprints and take this off-line or into chat to discuss. Please Google Scholar me and ping me. I also know the experimental folks and am happy to discuss. Jul 25 '21 at 2:32
• Let me put something coherent together. I'll also shoot a bunch of stuff up to github, which I should have done already. Wilson, Decius, and Cross is the canonical reference. There are matrix transformations (such as $\bf{R = LQ}$) where normal modes ($\bf{Q}$) are transformed to rectilinear combinations (local modes) which include $\ce{C-H}$ stretches, $\ce{H-C-C}$ bends (in-plane and out-of-plane), and so on. More later. Jul 25 '21 at 16:25
• For a much simpler intro, these systems can be treated semi-quantitatively as standard systems of coupled harmonic oscillators as summarized here: colby.edu/chemistry/PChem/notes/NormalModesText.pdf Jul 25 '21 at 16:45
• @ToddMinehardt Then the question title could use some improvement. I certainly would not expect a discussion of localized modes given the current title; also considering the lack of the term "localized modes" in the question body. It certainly would be interesting to see an explanation about localized modes though. I am also aware that there are a lot more details that could be discussed about normal modes, which are not given in the answer that I linked, but covering all details goes beyond the scope of one question, which is why I asked for clarification. Jul 26 '21 at 13:55
• @M.Farooq The answer to the old question does actually explain how you would obtain for example the ring breathing mode in benzene. I guess it isn't clear since that isn't the focus of that question/answer, but if you applied the described procedure to benzene, you would obtain the ring breathing mode. One of the $Q_i$ vectors would describe the motion that you have drawn in your figure. Localizing modes does something different. It's purpose is to convert a set of delocalized modes, like the ring breathing mode for example, to a set of localized modes that aren't spread out over the molecule Jul 26 '21 at 16:01

Vibrational spectroscopy and normal modes are based on the Born-Oppenheimer approximation and the harmonic approximation.

The BO approximation yields the adiabatic electronic states and the corresponding adiabatic energy eigenvalues as functions of the nuclear coordinates. The adiabatic energy eigenvalues of the electronic states are also known as potential energy surfaces $$E_\alpha(R)$$. Associated with this electronic adiabatic state is a vibronic Hamiltonian for the nuclear wave function, $$H_\alpha(R) = \frac{1}{2}\nabla^T_RM^{-1}\nabla_R + E_\alpha(R)$$

From now on I will fix the electronic state to be the groundstate and drop the index $$\alpha=0$$, that defines the electronic state. In pure vibrational spectroscopy the electronic state remains the same and is thus of no further interest.

If we knew the exact potential energy surface $$E(R)$$, then we could try to solve the problem numerically without the harmonic approximation. But doing so is very expensive for a large amount of degrees of freedom.

To avoid this, we can make the harmonic approximation and expand the potential energy surface to second order around a equilibrium structure, i.e. the nuclear coordinates where the potential energy surface has a minimum. Let this structure be $$R_{EQ}$$,

$$E(R)\approx E(R_{EQ}) + \sum^{3N}_{i=1}\frac{\partial E} {\partial R_i}|_{R_{EQ}} (R-R_{EQ})_i +\frac{1}{2}\sum^{3N,3N}_{i=1,j=1}(R-R_{EQ})_i \frac{\partial^2 E} {\partial^2 R_i R_j}|_{R_{EQ}} (R-R_{EQ})_j$$ We can chose our energy axis such that $$E(R_{EQ})=0$$. We also are at a minimum of the potential energy surface, which means that the gradient vanishes $$\frac{\partial E} {\partial R_i}|_{R_{EQ}}=0 \ \forall \ i$$. All that remains in the harmonic approximation at the equilbrium geometry is the second order term with the Hessian evaluated at the equilibrium geometry $$H_{ij}=\frac{\partial^2 E_\alpha} {\partial^2 R_i R_j}|_{R_{EQ}}$$, $$E^{Harmonic}(R) = \frac{1}{2}\sum^{3N,3N}_{i=1,j=1}(R-R_{EQ})_i H_{ij} (R-R_{EQ})_j$$

With this we can define the vibronic Hamiltonian in the harmonic approximation, \begin{aligned} H^{Harmonic}(R) &= \frac{1}{2}\nabla^T_RM^{-1}\nabla_R + E^{Harmonic}_\alpha(R)\\ H^{Harmonic}(R)&= \frac{1}{2}\nabla^T_RM^{-1}\nabla_R + \frac{1}{2}\sum^{3N,3N}_{i=1,j=1}(R-R_{EQ})_i H_{ij} (R-R_{EQ})_j \end{aligned}

A short remark on what I have not done is also in order. I have not separated the nuclear coordinates into translation, rotation and vibration thus far. Six(five for linear molecules) of the normal coordinates that we will obtain are going to describe translation and rotation. They are typically easy to recognize by their associated angular frequencies that are almost zero.

Now back to the problem. I will now drop the superscript harmonic. Our harmonic vibronic Hamilton operator is not a sum of 1D Harmonic oscillator Hamilton operators, but we can achieve that by doing some coordinate transformations.

First we introduce mass-weighted displacements $$x_i=\sqrt m_i (R-R_{EQ})_i \rightarrow \frac{\partial}{\partial x_i} = \frac{1}{\sqrt m_i}\frac{\partial}{\partial R_i}\rightarrow \nabla_x = M^{-1/2}\nabla_R$$ In these coordinates we obtain $$H(x) = \frac{1}{2}\nabla_x^T\nabla_x + \frac{1}{2}\sum^{3N,3N}_{i=1,j=1}x_i \frac{H_{ij}}{\sqrt{m_i m_j}} x_j$$ We define now the mass weighted Hessian $$F_{ij} = \frac{H_{ij}}{\sqrt{m_i m_j}}$$ which leads us to $$H(x) = \frac{1}{2}\nabla_x^T\nabla_x + \frac{1}{2}\sum^{3N,3N}_{i=1,j=1}x_i F_{ij} x_j$$

In this equation we still have coupling between different $$x_{i/j}$$ due to off-diagonal matrix elements in the mass weighted Hessian. But we can eliminate these by transformation to a coordinate system in which the matrix is diagonal. To find these coordinates we solve the eigenvalue problem $$FV = V\Lambda$$ where the matrix $$V$$ contains the eigenvectors of $$F$$ as columns. $$\Lambda$$ is the corresponding diagonal eigenvalue matrix.

The eigenvector matrix allows us to define the transformation to normal coordinates $$Q_i$$, $$\sum^{3N}_{k=1} V_{ik}Q_k = x_i \quad \rightarrow VQ=x \quad \leftrightarrow Q = V^{-1}x$$

The eigenvectors of the mass weighted Hessematrix are orthogonal and if normalized form a orthonormal set. This means that $$VV^T=1=VV^{-1}$$. Which means that $$F$$ can be diagonalized via $$\Lambda = V^TFV$$

Plugging in the transformation into the Hamiltonian yields \begin{aligned} H(Q) &= \nabla_Q^T({V^{-1}})^TV^{-1}\nabla_Q + \frac{1}{2} Q^TV^TFVQ \\ H(Q) &= \nabla_Q^T\nabla_Q +\frac{1}{2} Q^T \Lambda Q \end{aligned}

Expanding this equation shows that we now have a simple sum of harmonic 1D Hamiltonians since the Eigenvalues matrix is diagonal $$\Lambda_{ij}=\delta_{ij}\Lambda_i$$ and no longer couples different coordinates.The last step is recognizing that each 1D equation is the same as a harmonic oscillator equation with $$\omega_i^2=\Lambda_i.$$

\begin{aligned} H(Q) &= \sum^{3N}_{i=1}\frac{1}{2}\frac{\partial^2 }{\partial Q^2_i} + \sum^{3N}_{i=1}\frac{1}{2}Q_i^2\Lambda_i\\ &=\sum_{i=1}^{3N}\frac{1}{2}\frac{\partial^2}{\partial Q^2_i} + \frac{1}{2}Q_i^2\omega_i^2\\ H(Q) &= \sum_{i=1}^{3N} H^{1D}(Q_i) \end{aligned}

Meaning that our large coupled problem finally has taken the shape of a sum of 3N one dimensional harmonic oscillators, which we know how to solve analytically.

At this point I'd like to repeat that we have not separated translation and rotation. Due to this, the first six frequencies should be zero (or almost zero). The six normal modes corresponding to these zero eigenvalues are not oscillations and have to be omitted when calculating vibrational transitions.

It is possible to remove them before transforming to normal coordinates. How this is done is described in the book of Wilson, Decius and Cross, Molecular vibrations but discussing these points goes beyond the scope of my answer.

Let me now address how all of this connects to the figure that you show in your question. In particular how we know how the atoms are displaced relative to each other for a given normal mode.

By inverting the transformations, we see $$\Delta R = M^{-1/2}VQ$$ with $$\Delta R = R-R_{EQ}$$.

This means that the columns of the matrix \begin{aligned} K &= M^{-1/2}V \end{aligned} describe the relative cartesian displacement amplitudes of the normal modes, easily seen by plugging in a unit vector for Q.

If we had calculated benzene then one of the columns of $$K$$ would contain the relative displacements as shown in your figure. Some quantum chemistry programs also print these column vectors after some sort of normalization, when they output Cartesian normal modes and use them for visualization.

Now let me say a bit about reduced masses and force constants. First of all reduced mass and force constants are no observables. They are helpful quantities that we can define to make the analogy to a 1D system with a mass attached via spring to a wall. Depending on the particular conventions that we use, we will obtain different values. A common convention is to define a reduced mass based on the Cartesian displacement vector norm. The normal mode vectors as I have defined above are assumed to be normalized to one with regard to mass weighted coordinates. A mass weighted normal coordinate vector of norm one corresponds to a displacement of atoms such that all displacements squared, multiplied with their atomic mass then summed and the root taken yield 1. Take as an example one mode with many heavy atoms contributing to the motion and a mode where mostly light atoms contribute to motion. If both modes have now a normal coordinate value of 1, you would observe spatially when looking at the atom displacements, that the mode with the heavy atoms is only a little displaced while the mode with light atoms shows larger displacement.

We can however introduce a new coordinate that is constructed in a way that a value of 1 means always the same net spatial displacement of the atoms contributing to the mode irregardless of the mass of the contributing atoms.

Let us look at the norm of the Cartesian displacement vector associated with a mass weighted normal mode vector. A unit normal mode vector is of the form $$(\hat Q_\alpha)_j = \delta_{j\alpha}$$ The transformation of this vector to Cartesian displacement coordinates yields $$R_\alpha = K\hat Q_\alpha$$ The Cartesian norm squared that describes the actual displacement of the atoms in space is then, $$R_\alpha \cdot R_\alpha = K \hat Q_\alpha \cdot K\hat Q_\alpha = \sum_{nml} K_{nm} (\hat Q_\alpha)_m K_{nl} (\hat Q_\alpha)_l \\ = \sum_{nml} K_{nm} \delta_{\alpha m} K_{nl} \delta_{\alpha l} \\ = \sum_{nml} K_{n\alpha} K_{n\alpha} \\ |R_\alpha|^2 = \sum_{nml} (K_{n\alpha})^2 \\$$ We define now the reduced mass for mode $$\alpha$$ as $$\mu_\alpha = \frac{1}{|R_\alpha|^2}$$ and the new set of coordinates $$y_\alpha = \frac{1}{\sqrt \mu_\alpha}Q_\alpha$$ In this coordinates we obtain the Hamiltonians $$H(y_\alpha) = \sum_\alpha \frac{1}{2\mu_\alpha}\frac{\partial^2 }{\partial y_\alpha^2} + \frac{1}{2}\mu_\alpha \omega^2_\alpha y_\alpha^2$$ and can define the force constants for each mode $$k_\alpha = \mu_\alpha \omega^2_\alpha.$$

I am afraid you need a computer to do all of this or at least a very long breath and patience. I think for some systems the diagonalization can be done manually if you use symmetry, breaking down the matrices that have to be diagonalized to managable size, but generally you are going to let a computer do the linear algebra parts. I am pretty sure that is also the reason why almost all texts only discuss the trivial linear case where everything can be done by hand, although that is problematic since the simple case doesn't give you much intuition for the more general case.

All of this is also explained in the white paper https://gaussian.com/wp-content/uploads/dl/vib.pdf

I would not recommend the book from Wilson, Delcius and Cross as first reference. The book is good if you already have some understanding of the matter but as first exposition it might be difficult. It is however a good source if your are interested in internal coordinates and a treatment based on internal coordinates and the problem of separating vibrational and rotational degrees of freedom.