Classical Theory of Raman Scattering

When a system of molecules is placed in an electric field, the electrons and nuclei will be displaced in such a manner as to induce dipole moments. If we make the valid assumption that the electric field is weak, the induced dipole moment will be proportional to the field strength. The proportionality coefficient is called the polarizability, $\boldsymbol\alpha$ . Then in 3D space, the induced dipole moment is given by

$\boldsymbol{\mu} = \boldsymbol \alpha \cdot \boldsymbol{E}$

which is the dot product of a $3\times 3$ tensor and a $3$ dimentional column vector. If $\boldsymbol{E} = \boldsymbol{E}_0 \cos \omega t$ , then $\boldsymbol{\mu}$ will vary with time according to

$\boldsymbol\mu = \boldsymbol \alpha \cdot \boldsymbol E_0 \cos \omega t$

However, in real situations the vibration and the rotation of the molecule should be taken into account, which means that the polarizability is a function of time. For simplicity, let these be diatomic molecules and take small oscillation approximation:

$\boldsymbol \alpha = \boldsymbol \alpha_0 + \boldsymbol{\alpha_1}\cos \omega_0 t$

in which $\omega_0$ denotes for the frequency of the normal mode and $\alpha_1$ is a measure of how the polarizability varies with vibration.

Substitude, we have

$\boldsymbol \mu = \boldsymbol{\alpha_0}\cdot \boldsymbol{E_0} \cos \omega t + \boldsymbol{\alpha_1} \cdot \boldsymbol{E_0} \cos \omega t \cos \omega_0 t$

which further expands to

$\boxed{\boldsymbol{\mu} = \boldsymbol{\alpha_0} \cdot \boldsymbol{E_0} \cos \omega t + \frac{1}{2} \boldsymbol{\alpha_1}\cdot \boldsymbol{E_0} \cos (\omega + \omega_0) t + \frac{1}{2} \boldsymbol{\alpha_1}\cdot \boldsymbol{E_0} \cos (\omega - \omega_0) t}$

This equation shows us that the induced dipole oscillates not only with the incedent frequency $\omega$ , but also with the frequency $\omega \pm \omega_0$ . The first term accounts the Rayleigh scattering, and the second and third terms are called the Stokes and anti-Stokes scattering, respectively. We see that the frequencies observed in Raman scattering are beat frequencies between the radiation frequency $\omega$ and the molecular frequency $\omega_0$ , and the frequency difference between the incident light and the scattered light is called the Raman shift.

Classical Electrodynamics tells us that the mean rate of radiation emitted over all directions by a dipole oscillating according to $\boldsymbol{\mu} = \boldsymbol{\mu_0} \cos \omega t$ is given by

$I = \frac{\omega^4 }{3 c^3} \left| \boldsymbol{\mu_0} \right|^2$

where $c$ is the speed of light. The intensity of the scattered light is proportional to the square of the amplitude of the induced dipole moment, so there exist the proportionality where

$I_{\text{Raman}} \propto \left| \boldsymbol{\alpha_1} \cdot \boldsymbol{E_0} \right|^2 = \boldsymbol{E_0}^T \cdot \boldsymbol{\alpha_1}^T \cdot \boldsymbol{\alpha_1} \cdot \boldsymbol{E_0}$

and if some proper approximation is made (as we shall do later in quantum description of Raman scattering), the polarizability tensor could be symetric, thus we have

$I_{\text{Raman}} \propto \boldsymbol{E_0}^T \cdot \boldsymbol{\alpha_1}^2 \cdot \boldsymbol{E_0}$

Notes that $\boldsymbol{\alpha_1}$ stands for the small change of polarizability due to the vibration, so one can easily understand that when doing static calculation people always use the value $\left|\partial \boldsymbol{\alpha}/ \partial Q\right|^2$ at the equilibrium geometry, where $Q$ is the normal mode coordinate, to obtain the full spectroscopy. This was fully discussed in the 5th chapter of THE RAMAN EFFECT, 2002 by Derek A. Long.

Classical theory is not enough for us to understand the process of photon scattering by non relativistic atomic electrons. Concepts discussed above is just a rough physical picture of the Raman scattering process.

Quantum Theory of Raman Scattering

Kramers-Heisenberg Formula

With the help of Quantum Field Theory, one can treat the scattering process as creation and annihilation of photons and consider photons as quantum-mechanical excitations of the radiation field. The interaction Hamiltonian between the atomic electrons and the radation field is given by

$H_{\mathrm{int}}=\sum_i\left[-\frac{e}{2 m c}\left(\mathbf{p}_i \cdot \mathbf{A}\left(\mathbf{x}_i, t\right)+\mathbf{A}\left(\mathbf{x}_i, t\right) \cdot \mathbf{p}_i\right)+\frac{e^2}{2 m c^2} \mathbf{A}\left(\mathbf{x}_i, t\right) \cdot \mathbf{A}\left(\mathbf{x}_i, t\right)\right]$

where the summation is over the various atomic electrons that participate in the interaction, and the expression $\mathbf{A}\left(\mathbf{x}_i, t\right)$ originates from the vector potential and is now a field operator assumed to act on a photon state or a many-body state at $\mathbf{x}_i$ , where $\mathbf{x}_i$ is the position operator of the $i$ th electron.

Recalling the Time-Dependent Second-Order Perturbation Theory, the transition probability from an initial state $\left| i \right\rangle$ to a final state $\left| f \right\rangle$ is given by

$w_{i \rightarrow f}=\frac{2 \pi}{\hbar}\left|\left\langle f\left|\hat{H}_2\right| i\right\rangle+\sum_m \frac{\left\langle f\left|\hat{H}_1\right| m\right\rangle\left\langle m\left|\hat{H}_1\right| i\right\rangle}{E_f-E_i}\right|^2 \delta\left(E_f-E_i\right)$

If we describe the incedent light as with $(\hbar \omega, \boldsymbol\epsilon)$ where $\boldsymbol\epsilon$ is the unit vector of the incoming direction and also the scattered light $(\hbar\omega', \boldsymbol\epsilon')$ where $\boldsymbol\epsilon'$ is the unit vector of outcoming direction, after some fancy derivation in J. J. Sakurai’s book ADVANCED QUANTUM MECHANICS, 1967, the transition probability in some solid angle can be written as

$\begin{gathered} w_{d \Omega}=\frac{2 \pi}{\hbar}\left(\frac{c^2 \hbar}{2 V \sqrt{\omega \omega^{\prime}}}\right)^2\left(\frac{e^2}{m c^2}\right)^2 \frac{V}{(2 \pi)^3} \frac{\omega^{\prime 2}}{\hbar c^3} d \Omega \\ \times\left|\delta_{i f} \boldsymbol{\epsilon} \cdot \boldsymbol{\epsilon}'-\frac{1}{m} \sum_n\left(\frac{\langle f \left|\mathbf{\boldsymbol{\hat{\mu}}} \cdot \boldsymbol{\epsilon}'\right|n\rangle \langle n \left| \mathbf{\boldsymbol{\hat{\mu}}} \cdot \boldsymbol{\epsilon}\right|i\rangle}{E_n-E_i-\hbar \omega}-\frac{\langle f \left|\mathbf{\boldsymbol{\hat{\mu}}} \cdot \boldsymbol{\epsilon} \right|n \rangle \langle n \left| \mathbf{\boldsymbol{\hat{\mu}}} \cdot \boldsymbol{\epsilon}'\right| i \rangle }{E_n-E_i+\hbar \omega^{\prime}}\right)\right|^2\delta\left(E_f-E_i\right) \end{gathered}$

where $\boldsymbol{\hat{\mu}}$ is the dipole operator (sum of charge times position) and $r_0$ stands for the classical radius of the electron

$r_0=\frac{e^2}{4 \pi m c^2} \simeq \frac{1}{137} \frac{\hbar}{m c} \simeq 2.82 \times 10^{-13} \mathrm{~cm}$

When divided by the flux density $c/V$ and solid angle $d\Omega$ , the transition probability becomes the differential cross section:

$\frac{d \sigma}{d \Omega}=r_0^2\left(\frac{\omega^{\prime}}{\omega}\right) \left|\delta_{i f} \boldsymbol{\epsilon} \cdot \boldsymbol{\epsilon}'-\frac{1}{m} \sum_n\left(\frac{\langle f \left|\mathbf{\boldsymbol{\hat{\mu}}} \cdot \boldsymbol{\epsilon}'\right|n\rangle \langle n \left| \mathbf{\boldsymbol{\hat{\mu}}} \cdot \boldsymbol{\epsilon}\right|i\rangle}{E_n-E_i-\hbar \omega}-\frac{\langle f \left|\mathbf{\boldsymbol{\hat{\mu}}} \cdot \boldsymbol{\epsilon} \right|n \rangle \langle n \left| \mathbf{\boldsymbol{\hat{\mu}}} \cdot \boldsymbol{\epsilon}'\right| i \rangle }{E_n-E_i+\hbar \omega^{\prime}}\right)\right|^2\delta\left(E_f-E_i\right)$

This is called the Kramers-Heisenberg formula. It was derived before the advent of quantum mechanics by Hendrik Kramers and Werner Heisenberg in 1925, based on the correspondence principle applied to the classical dispersion formula for light. The quantum mechanical derivation was given by Paul Dirac in 1927. When the final state and the initial state do not coincide (as always the case in Raman scattering), the first term in the square bracket vanishes, and because the scattered frequency is usually close to the incident frequency, $\omega \simeq \omega'$ , we may write the differential corss section of Raman scattreing as

$\frac{d \sigma}{d \Omega}= \frac{r_0^2}{m} \left|\langle i \left| \boldsymbol{\epsilon} \cdot \boldsymbol{\hat{\alpha }} \cdot \boldsymbol{\epsilon}'\right| f \rangle \right|^2\delta\left(E_f-E_i\right),\quad \text{where}\; \boldsymbol{\hat{\alpha}} = \sum_n\left[\frac{\mathbf{\boldsymbol{\hat{\mu}}}\, |n\rangle \langle n | \,\mathbf{\boldsymbol{\hat{\mu}}}}{E_n-E_i-\hbar \omega}-\frac{\mathbf{\boldsymbol{\hat{\mu}}} \,|n \rangle \langle n |\, \mathbf{\boldsymbol{\hat{\mu}}} }{E_n-E_i+\hbar \omega}\right]$

differentiate $\displaystyle\frac{d \sigma}{d \Omega}$ with respect to $\hbar \omega$ , we have

$\boxed{\frac{d }{d\hbar\omega} \frac{d \sigma}{d \Omega}= \frac{r_0^2}{m} \left|\langle i \left| \boldsymbol{\epsilon} \cdot \boldsymbol{\hat{\alpha }} \cdot \boldsymbol{\epsilon}'\right| f \rangle \right|^2\delta\left(\omega - \omega_{fi}\right), \quad \text{where}\; \omega_{fi} = \frac{E_f - E_i}{\hbar}}$

This could be the starting point of introducing the correlation function formalism for AIMD simulation, but there exists one main obstacle: the dependence on $n$ in the denominator of $\boldsymbol{\hat{\alpha}}$ . So we have to do some simplification to make the formula more tractable.

Simplification for the Polarizability

Learned from the book RAMAN SPECTROSCOPY: THEORY AND PRACTICE by Herman A. Szymanski, this siplification is straight forward, involving a series expansion and a truncation. Let’s take out the first term

$\boldsymbol{\hat{\alpha}_1} = \sum_n\frac{\mathbf{\boldsymbol{\hat{\mu}}}\, |n\rangle \langle n | \,\mathbf{\boldsymbol{\hat{\mu}}}}{E_n-E_i-\hbar \omega}$

as an example and later we will find that this mathematical trick suits the second term $\boldsymbol{\hat{\alpha}_2}$ as well.

Set the initial state always be the ground state, so we substitude $E_0$ for $E_i$ . I claim that there EXIST some constant "average energy" $E_{av}$ such that the following equation holds:

$\boxed{\frac{1}{E_n - E_0 - \hbar \omega} = \sum_{N=0}^\infty (-1)^N \frac{(E_n-E_{av})^N}{(E_{av} - E_0 - \hbar \omega)^{N+1}}, \quad \text{if}\; \left|\frac{E_n-E_{av}}{E_{av} - E_0 - \hbar \omega}\right| < 1}$

This is in fact a geometric series expansion, and the requirement behind is the convergence condition of the geometric series. Plug into $\boldsymbol{\hat{\alpha}_1}$ , we have

$\boldsymbol{\hat{\alpha}_1}= \sum_n\frac{\mathbf{\boldsymbol{\hat{\mu}}}\, |n\rangle \langle n | \,\mathbf{\boldsymbol{\hat{\mu}}}}{E_n-E_0-\hbar \omega} = \sum_n\sum_{N=0}^\infty (-1)^N \mathbf{\boldsymbol{\hat{\mu}}}\left[\frac{(E_n-E_{av})^N}{(E_{av} - E_0 - \hbar \omega)^{N+1}}\right] |n\rangle \langle n | \,\mathbf{\boldsymbol{\hat{\mu}}}$

Because $E_n$ is the eigenvalue of the Hamiltonian, we can substitude $\hat{H}|n\rangle$ for $E_n |n\rangle$ in the above equation:

$\boldsymbol{\hat{\alpha}_1}= \sum_n\sum_{N=0}^\infty (-1)^N \mathbf{\boldsymbol{\hat{\mu}}}\left[\frac{(\hat{H}-E_{av})^N}{(E_{av} - E_0 - \hbar \omega)^{N+1}}\right] |n\rangle \langle n | \,\mathbf{\boldsymbol{\hat{\mu}}}$

swich the order of the summation, we eliminate the summation over $n$ and obtain

$\boldsymbol{\hat{\alpha}_1}= \sum_{N=0}^\infty (-1)^N \frac{\mathbf{\boldsymbol{\hat{\mu}}}\,(\hat{H}-E_{av})^N\,\mathbf{\boldsymbol{\hat{\mu}}}}{(E_{av} - E_0 - \hbar \omega)^{N+1}} = \frac{\mathbf{\boldsymbol{\hat{\mu}}}\,\mathbf{\boldsymbol{\hat{\mu}}}}{E_{av} - E_0 - \hbar \omega} - \frac{\mathbf{\boldsymbol{\hat{\mu}}}\,(\hat{H}-E_{av})\,\mathbf{\boldsymbol{\hat{\mu}}}}{(E_{av} - E_0 - \hbar \omega)^2} + \cdots$

Note that this series is not only a geometric series but also alternating series, so it converges pretty fast. If we truncate the series at the first term when $N = 0$ and assume the polarizability does not change too much when light is shed on the molecule (which means this is the zero-order approximation), we have the following proportionality:

$\boldsymbol{\hat{\alpha}_1} \propto \mathbf{\boldsymbol{\hat{\mu}}}\,\mathbf{\boldsymbol{\hat{\mu}}} \quad \text{or} \quad \hat{\alpha}_{1,\rho\sigma} \propto \mathbf{\hat{\mu}}_\rho\,\mathbf{\hat{\mu}}_\sigma$

where $\rho,\sigma$ run from 1 to 3 denoting for $x,y,z$ in 3D space.

Obviously, this mathematical trick is not only applicable to the first term $\boldsymbol{\hat{\alpha}_1}$ , but also to the second term $\boldsymbol{\hat{\alpha}_2}$ . Put them together, what really matters is the following proportionality:

$\boxed{\boldsymbol{\hat{\alpha}} \propto \mathbf{\boldsymbol{\hat{\mu}}}\,\mathbf{\boldsymbol{\hat{\mu}}} \propto \mathbf{\boldsymbol{\hat{x}}}\,\mathbf{\boldsymbol{\hat{x}}} \quad \text{or} \quad \hat{\alpha}_{\rho\sigma} \propto \mathbf{\hat{\mu}}_\rho\,\mathbf{\hat{\mu}}_\sigma \propto \hat{x}_\rho\,\hat{x}_\sigma}$

By the elegance of proportionality, I have eliminated the somehow weird $E_{av}$ . The EXISTENCE of $E_{av}$ is all I need! Go back, we may now introduce our correlation function formalism.

Correlation Function formalism

This part follows Chapter 21: The Time-Correlation Function Formalism I, STATISTICAL MECHANICS, 1976 by Donald A. McQuarrie. Do ensemble average to the differential cross section and switch the delta function to the time domain, we have

$\frac{d }{d\hbar\omega} \frac{d \sigma}{d \Omega} \propto \int_{-\infty}^{\infty}dt e^{-i \omega t}\sum_{if}\rho_i \left|\langle i \left| \boldsymbol{\epsilon} \cdot \boldsymbol{\hat{\alpha }} \cdot \boldsymbol{\epsilon}'\right| f \rangle \right|^2 e^{i(\omega_f - \omega_i)t}$

where $\rho_i$ is the density distribution of the initial state. It is pretty easy to see that the summation over $i,f$ gives ensemble average, and handing out the $e^{i\omega_f t}, e^{ - i\omega_it}$ paves the way to Dirac interaction picture, which gives

$\frac{d }{d\hbar\omega} \frac{d \sigma}{d \Omega} \propto \int_{-\infty}^{\infty}dt e^{-i \omega t}\left\langle (\boldsymbol{\epsilon} \cdot \boldsymbol{\hat{\alpha}}(0) \cdot \boldsymbol{\epsilon}' )(\boldsymbol{\epsilon} \cdot \boldsymbol{\hat{\alpha}}(t) \cdot \boldsymbol{\epsilon}' )\right\rangle$

where $(\boldsymbol{\epsilon} \cdot \boldsymbol{\hat{\alpha}}(0) \cdot \boldsymbol{\epsilon}' )(\boldsymbol{\epsilon} \cdot \boldsymbol{\hat{\alpha}}(t) \cdot \boldsymbol{\epsilon}' )$ can be written in component form as

$(\boldsymbol{\epsilon} \cdot \boldsymbol{\hat{\alpha}}(0) \cdot \boldsymbol{\epsilon}' )(\boldsymbol{\epsilon} \cdot \boldsymbol{\hat{\alpha}}(t) \cdot \boldsymbol{\epsilon}' ) = \left(\sum_{i=1}^3\sum_{j=1}^3\epsilon_i\, \hat\alpha_{ij}(0)\,\epsilon_j'\right) \left(\sum_{k=1}^3\sum_{l=1}^3\epsilon_k\,\hat\alpha_{kl}(t)\,\epsilon_l'\right)$

in which $\alpha_{ij}(t)$ is the $ij$ th component of the polarizability tensor. Our differential cross section is now

$\frac{d }{d\hbar\omega} \frac{d \sigma}{d \Omega} \propto \int_{-\infty}^{\infty}dt e^{-i \omega t}\sum_{ijkl}\left\langle \epsilon_i\epsilon_j'\epsilon_k\epsilon_l'\right\rangle_{\text{sphere}} \times \left\langle \hat\alpha_{ij}(0)\hat\alpha_{kl}(t)\right\rangle_{\text{ensemble}}$

where $\left\langle \epsilon_i\epsilon_j'\epsilon_k\epsilon_l'\right\rangle_{\text{sphere}}$ is the average over the solid angle. In experiments, people always put the observer perpendicular to the incident light beam, so we set $\boldsymbol{\epsilon} \perp \boldsymbol{\epsilon}'$ , thus we have

$\epsilon_i\epsilon_j'\epsilon_k\epsilon_l' = \frac{1}{30} (4\delta_{ij}\delta_{kl} - \delta_{ik}\delta_{jl} - \delta_{il}\delta_{jk})$

Plug this, it becomes

$\frac{d }{d\hbar\omega} \frac{d \sigma}{d \Omega} \propto \int_{-\infty}^{\infty}dt e^{-i \omega t}\sum_{\rho\sigma}\left\langle \hat\alpha_{\rho\sigma}(0)\hat\alpha_{\sigma\rho}(t) - \frac{1}{3} \hat\alpha_{\rho\rho}(0)\hat\alpha_{\sigma\sigma}(t) \right\rangle_{\text{ensemble}}$

which can be written in a more compact form:

$\boxed{\frac{d }{d\hbar\omega} \frac{d \sigma}{d \Omega} \propto \int_{-\infty}^{\infty}dt e^{-i \omega t}\left\langle \text{Tr}\left[\boldsymbol{\hat{\beta}}(0) \cdot \boldsymbol{\hat{\beta}}(t)\right]\right\rangle_{\text{ensemble}}}$

where $\boldsymbol{\hat{\beta}}(t) = \boldsymbol{\hat{\alpha}}(t) - \frac{1}{3}\text{Tr}\left[\boldsymbol{\hat{\alpha}}(t)\right]\boldsymbol{I}$ , so that $\text{Tr}[\boldsymbol{\hat{\beta}}(t)] = 0$ .

With the help of our simplification in section 2.2, we can repalce the polarizability tensor with the position operator to give

$\frac{d }{d\hbar\omega} \frac{d \sigma}{d \Omega} \propto \int_{-\infty}^{\infty}dt e^{-i \omega t}\sum_{\rho\sigma}\left\langle \hat{x}_\rho(0)\hat{x}_\sigma(0)\hat{x}_\sigma(t)\hat{x}_\rho(t) - \frac{1}{3}\, \hat{x}_\rho(0)\hat{x}_\rho(0)\hat{x}_\sigma(t)\hat{x}_\sigma(t) \right\rangle_{\text{ensemble}}$

In practice, because quantum correlation functions are usually difficult to calculate, classical autocorrelation functions obtained from classical simulations are more often used for spectra computation. Within the harmonic oscillator model, the Fourier transform of quantum and classical autocorrelation function only differ by a frequency-dependent prefactor, as we have already seen in infrared spectroscopy.

Although the idea for calculating both classical and quantum average is quite simple, as is mentioned in Yiwen’s note for IR spectrum, summation over $\rho,\sigma$ in Raman spectrum brings more difficulties above. Introducing 3D isotropic harmonic approximation:

$H = \frac{p_1^2+ p_2^2 + p_3^2}{2m} + \frac{1}{2} m \omega_0^2 (x_1^2 + x_2^2 + x_3^2)$

thanks to Wolfram Mathematica 13.0.0, Sthdent Edition, the classical average takes

$\frac{10 k^2T^2}{m^2 \omega_0^4}\cos^2\omega_0 t$

and the quantum average takes

$\frac{\hbar ^2 e^{\frac{\left(3 k^2T^2-1\right) \hbar \omega_0 }{2 kT}} \left[7 \cos \left(2 w_0 t-\frac{i \hbar \omega_0 }{kT}\right)-3 \cosh \left(\frac{\hbar\omega_0 }{kT}\right)+3 \cos (2 \omega_0 t)+13\right]}{2 m^2 \omega_0^2 \left(e^{\frac{\hbar\omega_0}{kT}}-1\right)^2}$

After Fourier transformation, the prefactor is determined by

$\boxed{\int_{-\infty}^{\infty}dt e^{-i \omega t} \left\langle \cdots \right\rangle_{\text{Quantum}} \;\Big/\int_{-\infty}^{\infty}dt e^{-i \omega t} \left\langle \cdots \right\rangle_{\text{Classical}} = \frac{e^{(-1+3k^2T^2)\hbar \omega_0/2kT}(3+7e^{\hbar \omega_0 / kT})\hbar^2\omega_0^2}{10(-1+e^{\hbar \omega_0 / kT})^2 k^2T^2}}$

The Mathematica .nb file is attached as appendix. Please contact me if
you find any mistakes.

Conclusion

In this note, I start with the classical theory for light scattering, giving a physical picture of this process and then terned to quantum description. The Kramers-Heisenberg Formula comes from the Time-Dependent Second-Order Perturbation Theory which focuses on the transition probability of a system under harmonic perturbation. As the transition probability is proportional to the intensity of scattered light, we then jump to the expression of scattering differential cross section, where we begin the rest of our journey.

The correlation function formalism cannot be introduced directly because the summation over states in the oringial formula. Under zero-th-order approximation, we let some “average energy” of the states to replace the summation over states, and then siplification of the polarizability leads to the proportionality between the polarizability tensor and the summation and multiplication of position operators, which is further calculated by mathematical software.

One thing I have to emphasize is that the truncation of our polarizability and later the introduction of proportionality together get rid of the frequency dependence ( $\omega$ in the denominator) of the original polarizability. This only holds when the polarizability does not change too much during the measurement. \textcolor{RoyalBlue}{\textbf{In some situations, this could have some influence on the final prefactor}}. Let $f_{\text{q}}(\lambda)$ and $f_{\text{c}}(\lambda)$ and stands for the frequency dependence of the polarizability, then the prefactor could look like

$\boxed{\int_{-\infty}^{\infty}dt e^{-i \omega t} \left\langle \cdots \right\rangle_{\text{Quantum}} \;\Big/\int_{-\infty}^{\infty}dt e^{-i \omega t} \left\langle \cdots \right\rangle_{\text{Classical}} = \frac{\displaystyle\int f_{\text{q}}(\lambda) d \lambda }{\displaystyle\int f_{\text{c}}(\lambda) d \lambda } \frac{e^{(-1+3k^2T^2)\hbar \omega_0/2kT}(3+7e^{\hbar \omega_0 / kT})\hbar^2\omega_0^2}{10(-1+e^{\hbar \omega_0 / kT})^2 k^2T^2}}$