Skip to main content
Chemistry LibreTexts

5.I1: Simple Harmonic Oscillator - Plotting Eigenstates

  • Page ID
    281077
  • \( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

    \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)

    \( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)

    ( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)

    \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)

    \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)

    \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)

    \( \newcommand{\Span}{\mathrm{span}}\)

    \( \newcommand{\id}{\mathrm{id}}\)

    \( \newcommand{\Span}{\mathrm{span}}\)

    \( \newcommand{\kernel}{\mathrm{null}\,}\)

    \( \newcommand{\range}{\mathrm{range}\,}\)

    \( \newcommand{\RealPart}{\mathrm{Re}}\)

    \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)

    \( \newcommand{\Argument}{\mathrm{Arg}}\)

    \( \newcommand{\norm}[1]{\| #1 \|}\)

    \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)

    \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)

    \( \newcommand{\vectorA}[1]{\vec{#1}}      % arrow\)

    \( \newcommand{\vectorAt}[1]{\vec{\text{#1}}}      % arrow\)

    \( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

    \( \newcommand{\vectorC}[1]{\textbf{#1}} \)

    \( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)

    \( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)

    \( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)

    \( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

    \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)

    \(\newcommand{\avec}{\mathbf a}\) \(\newcommand{\bvec}{\mathbf b}\) \(\newcommand{\cvec}{\mathbf c}\) \(\newcommand{\dvec}{\mathbf d}\) \(\newcommand{\dtil}{\widetilde{\mathbf d}}\) \(\newcommand{\evec}{\mathbf e}\) \(\newcommand{\fvec}{\mathbf f}\) \(\newcommand{\nvec}{\mathbf n}\) \(\newcommand{\pvec}{\mathbf p}\) \(\newcommand{\qvec}{\mathbf q}\) \(\newcommand{\svec}{\mathbf s}\) \(\newcommand{\tvec}{\mathbf t}\) \(\newcommand{\uvec}{\mathbf u}\) \(\newcommand{\vvec}{\mathbf v}\) \(\newcommand{\wvec}{\mathbf w}\) \(\newcommand{\xvec}{\mathbf x}\) \(\newcommand{\yvec}{\mathbf y}\) \(\newcommand{\zvec}{\mathbf z}\) \(\newcommand{\rvec}{\mathbf r}\) \(\newcommand{\mvec}{\mathbf m}\) \(\newcommand{\zerovec}{\mathbf 0}\) \(\newcommand{\onevec}{\mathbf 1}\) \(\newcommand{\real}{\mathbb R}\) \(\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}\) \(\newcommand{\laspan}[1]{\text{Span}\{#1\}}\) \(\newcommand{\bcal}{\cal B}\) \(\newcommand{\ccal}{\cal C}\) \(\newcommand{\scal}{\cal S}\) \(\newcommand{\wcal}{\cal W}\) \(\newcommand{\ecal}{\cal E}\) \(\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}\) \(\newcommand{\gray}[1]{\color{gray}{#1}}\) \(\newcommand{\lgray}[1]{\color{lightgray}{#1}}\) \(\newcommand{\rank}{\operatorname{rank}}\) \(\newcommand{\row}{\text{Row}}\) \(\newcommand{\col}{\text{Col}}\) \(\renewcommand{\row}{\text{Row}}\) \(\newcommand{\nul}{\text{Nul}}\) \(\newcommand{\var}{\text{Var}}\) \(\newcommand{\corr}{\text{corr}}\) \(\newcommand{\len}[1]{\left|#1\right|}\) \(\newcommand{\bbar}{\overline{\bvec}}\) \(\newcommand{\bhat}{\widehat{\bvec}}\) \(\newcommand{\bperp}{\bvec^\perp}\) \(\newcommand{\xhat}{\widehat{\xvec}}\) \(\newcommand{\vhat}{\widehat{\vvec}}\) \(\newcommand{\uhat}{\widehat{\uvec}}\) \(\newcommand{\what}{\widehat{\wvec}}\) \(\newcommand{\Sighat}{\widehat{\Sigma}}\) \(\newcommand{\lt}{<}\) \(\newcommand{\gt}{>}\) \(\newcommand{\amp}{&}\) \(\definecolor{fillinmathshade}{gray}{0.9}\)

    The Harmonic Oscillator (HO) is one of the most important systems in quantum mechanics for the following reasons:

    1. It can solved analytically
    2. A harmonic potential can be used to approximate many other, more complex potentials
    3. It is the foundation for more advanced physics, from understanding phonons in materials to quantum field theory

    Question 1: Hermite Polynomials

    We know that the time independent, normalized stationary states for the HO are given in terms of Hermite polynomials, \(H_{n}(\xi)\), by

    \[\xi=\sqrt{\frac{m \omega}{\hbar}} x\]

    import matplotlib
    import matplotlib.pyplot as plt 
    import numpy
    import numpy.polynomial.hermite as Herm
    
    #Choose simple units
    m=1.
    w=1.
    hbar=1.
    
    #Discretized space
    dx = 0.05
    x_lim = 12
    x = numpy.arange(-x_lim,x_lim,dx)
    
    def hermite(x, n):
        xi = numpy.sqrt(m*w/hbar)*x
        herm_coeffs = numpy.zeros(n+1)
        herm_coeffs[n] = 1
        return Herm.hermval(xi, herm_coeffs)
      
    plt.figure()
    plt.plot(x, hermite(x,0), linewidth=2,label=r"$H_0$")
    plt.plot(x, hermite(x,1), linewidth=2,label=r"$H_1$")
    plt.plot(x, hermite(x,2), linewidth=2,label=r"$H_2$")
    plt.plot(x, hermite(x,3), linewidth=2,label=r"$H_3$")
    plt.plot(x, hermite(x,4), linewidth=2,label=r"$H_4$")
    
    #Set limits for axes
    plt.xlim([-2.5,2.5])
    plt.ylim([-20,20])
    
    #Set axes labels
    plt.xlabel("x")
    plt.ylabel(r"$H_n(\xi)$")
    plt.title(r"Hermite Polynomials, $H_n(\xi)$")
    plt.legend()
    plt.show()

    The Harmonic Oscillator Eigenstates

    import matplotlib
    import matplotlib.pyplot as plt 
    import numpy
    import numpy.polynomial.hermite as Herm
    import math
    
    #Choose simple units
    m=1.
    w=1.
    hbar=1.
    #Discretized space
    dx = 0.05
    x_lim = 12
    x = numpy.arange(-x_lim,x_lim,dx)
    
    def hermite(x, n):
        xi = numpy.sqrt(m*w/hbar)*x
        herm_coeffs = numpy.zeros(n+1)
        herm_coeffs[n] = 1
        return Herm.hermval(xi, herm_coeffs)
      
    def stationary_state(x,n):
        xi = numpy.sqrt(m*w/hbar)*x
        prefactor = 1./math.sqrt(2.**n * math.factorial(n)) * (m*w/(numpy.pi*hbar))**(0.25)
        psi = prefactor * numpy.exp(- xi**2 / 2) * hermite(x,n)
        return psi
    
    plt.figure()
    plt.plot(x, stationary_state(x,4))
    plt.xlabel(r"x")
    plt.ylabel(r"$\psi_4(x)$")
    plt.title(r"Test Plot of $\psi_4(x)$")
    plt.show()

    The Classical Harmonic Oscillator

    Now returns the classical probability density for a particle with an energy equal to that of the $n^{th}$ quantum harmonic oscillator energy level.

    import matplotlib
    import matplotlib.pyplot as plt 
    import numpy
    import numpy.polynomial.hermite as Herm
    
    #Choose simple units
    m=1.
    w=1.
    hbar=1.
    #Discretized space
    dx = 0.05
    x_lim = 12
    x = numpy.arange(-x_lim,x_lim,dx)
    
    def classical_P(x,n):
        E = hbar*w*(n+0.5)
        x_max = numpy.sqrt(2*E/(m*w**2))
        classical_prob = numpy.zeros(x.shape[0])
        x_inside = abs(x) < (x_max - 0.025)
        classical_prob[x_inside] = 1./numpy.pi/numpy.sqrt(x_max**2-x[x_inside]*x[x_inside])
        return classical_prob
    
    plt.figure()
    plt.plot(x, classical_P(x,4))
    plt.xlabel(r"x")
    plt.ylabel(r"$P_{classical}(x,4)$")
    plt.title(r"Classical Probability Density of a Harmonic Oscillator, $E=E_4$")
    plt.show()

    Comparing Classical vs. Quantum Harmonic Results

    Compare the behavior of a quantum harmonic oscillator to a classical harmonic oscillator. Connect what happens as you increase the quantum number to the transition from quantum to classical behavior. There it is shown that for a classical harmonic oscillator with energy \(E\), the classical probability of finding the particle at \(x\) is given by

    \[x_{\max }=\sqrt{\frac{2 E}{m \omega^{2}}}\]

    where \(x_{\max }\) is the classical turning point, and \(P_{\text {classical}}(x)\) is understood to be zero for \(|x|>\left|x_{\max }\right|\).

    import matplotlib
    import matplotlib.pyplot as plt 
    import numpy
    import numpy.polynomial.hermite as Herm
    import math
    
    #Choose simple units
    m=1.
    w=1.
    hbar=1.
    #Discretized space
    dx = 0.05
    x_lim = 12
    x = numpy.arange(-x_lim,x_lim,dx)
    
    def hermite(x, n):
        xi = numpy.sqrt(m*w/hbar)*x
        herm_coeffs = numpy.zeros(n+1)
        herm_coeffs[n] = 1
        return Herm.hermval(xi, herm_coeffs)
    def stationary_state(x,n):
        xi = numpy.sqrt(m*w/hbar)*x
        prefactor = 1./math.sqrt(2.**n * math.factorial(n)) * (m*w/(numpy.pi*hbar))**(0.25)
        psi = prefactor * numpy.exp(- xi**2 / 2) * hermite(x,n)
        return psi
      
    def classical_P(x,n):
        E = hbar*w*(n+0.5)
        x_max = numpy.sqrt(2*E/(m*w**2))
        classical_prob = numpy.zeros(x.shape[0])
        x_inside = abs(x) < (x_max - 0.025)
        classical_prob[x_inside] = 1./numpy.pi/numpy.sqrt(x_max**2-x[x_inside]*x[x_inside])
        return classical_prob
    
    plt.figure(figsize=(10, 8))
    plt.subplot(3,2,1)
    plt.plot(x, numpy.conjugate(stationary_state(x,0))*stationary_state(x,0), label="n=0")
    plt.plot(x, classical_P(x,0))
    plt.legend()
    plt.subplot(3,2,2)
    plt.plot(x, numpy.conjugate(stationary_state(x,3))*stationary_state(x,3), label="n=3")
    plt.plot(x, classical_P(x,3))
    plt.legend()
    plt.subplot(3,2,3)
    plt.plot(x, numpy.conjugate(stationary_state(x,8))*stationary_state(x,8), label="n=8")
    plt.plot(x, classical_P(x,8))
    plt.legend()
    plt.subplot(3,2,4)
    plt.plot(x, numpy.conjugate(stationary_state(x,15))*stationary_state(x,15), label="n=15")
    plt.plot(x, classical_P(x,15))
    plt.legend()
    plt.subplot(3,2,5)
    plt.plot(x, numpy.conjugate(stationary_state(x,25))*stationary_state(x,25), label="n=25")
    plt.plot(x, classical_P(x,25))
    plt.legend()
    plt.subplot(3,2,6)
    plt.plot(x, numpy.conjugate(stationary_state(x,40))*stationary_state(x,40), label="n=40")
    plt.plot(x, classical_P(x,40))
    plt.legend()
    plt.show()

    Q5: Comparison of the Classical and Quantum Harmonic Oscillators

    We know that quantum physics often makes counter-intuitive predictions that contradict classical physics. But we also know that quantum effects do not often appear in the macroscopic world. For the example of the harmonic oscillator, at which energies do classical and quantum physics most agree, and at which energies do they most strongly diverge? Based on the above graphs, describe some of the differences predicted by quantum and classical physics for the harmonic oscillator.

    The difference between quantum and classical physics is most distinct for the lower energy levels, where quantum mechanics does not at all obey the classical intuition that particles spend most of their time near the classical turning point. At very high energies though, the general shape of the quantum probability distribution starts to look more like the classical result.

    Here are three additional differences between the classical and quantum results.

    1. Classically, there is zero probability for the particle to exist outside of the energetically allowed region (beyond the turning point). However, the quantum results indicate some small but finite probability of measuring the particle in the classically "forbidden" region.
    2. The classical probability distribution says that there is some chance of finding the particle at any point within the classically allowed region. In contrast, the quantum probability distributions possess nodes of zero probability within this region.
    3. This one's not obvious from the graphs, but it's clear enough that a classical harmonic oscillator can possess any energy. In particular, there is no minimum allowable energy, in stark contrast to the quantum harmonic oscillator, whose minimum energy (ground state energy, vacuum energy) is \(E_{0}=\hbar \omega / 2\).

    5.I1: Simple Harmonic Oscillator - Plotting Eigenstates is shared under a not declared license and was authored, remixed, and/or curated by LibreTexts.

    • Was this article helpful?