Quantum Harmonic Oscillator.sagews

Physicists LOVE the harmonic oscillator potential. It is one of our favorites for two very important reasons: 1) It is one of handful of potentials that we can solve exactly and 2) any restoring force potential can be approximated as a harmonic oscillator near the equilibrium point. This is why we often approximate solid matter as a collection of ball-and-springs; the molecular bonds between atoms behave like springs if you don't compress the solid too much.

The harmonic potential energy equation is (U=\frac{1}{2}m \omega^2 x^2). Plugging this into Schrodinger's equation gives

[\frac{-\hbar^2}{2m} \frac{d^2}{dx^2}\psi(x) + \frac{1}{2}m \omega^2 x^2 \psi(x) = E \psi(x)]

Solving differential equations

• When solving differential equations it helps to specify which variables are real and which are positive using the assume() command.

• You need to specify the independent variable of your function. In this case we specify that $\psi$ is a function of $x$.

• The differential equation is defined using the Python def method.

  • This creates a function which we can differential, integrate, or solve using the differential equation solver.

• Use desolve() to solve the differential equation. The arguments of desolve() are the differential equation to solve, the dependent variable ($\psi$ in this case), and the independent variable x (specified using ivar=x to indicate the independent variable). In this case, desolve() with encounter an issue unless you specify contrib_ode=True, in which case a slightly different solver is used to handle the differential equation.

psi,m,L,omega=var('psi,m,L,omega')  #Declare wave function, mass, and box length as variables

Ener=var('Ener',latex_name=r'\mathcal{E}')  #Declare Ener as energy variable but will display a script E when using show() and pretty_print
hbar=var('hbar',latex_name=r'\hbar')  #Declares hbar as a variable and displays hbar symbol when using show() and pretty_print

assume(Ener,'real')   #Assume the variable Ener (the energy) only takes real values...
assume(Ener>0)   #... and assume the energy is greater than zero.  These are needed for desolve()

assume(m,'real')   #Assume the mass is real...
assume(m>0)     #... and positive.  If you don't make these assumptions you will get an error message asking if these variables are positive, negative, or zero

psi=function('psi')(x)   #To solve the differential equation we must tell SageMath that psi is a function called 'psi' with a single dependent variable


def f(psi):  #Define the equation to solve.
    return -hbar^2/(2*m)*diff(psi,x,2)+(1/2*m*omega^2*x^2-Ener)*psi   #This is the differential equation to be solved

general_soln = desolve(f(psi)==0,psi,ivar=x,contrib_ode=True)   #


pretty_print('The general solution to Schrodinger\'s equation is ', general_soln[0])

The general solution to Schrodinger's equation is $\displaystyle \psi\left(x\right) = K_{1} e^{\left(\frac{m \omega x^{2}}{2 \, {{\hbar}}}\right)} {\rm kummer}_{m}\left(\frac{{{\hbar}} \omega + 2 \, {{\mathcal{E}}}}{4 \, {{\hbar}} \omega}, \frac{1}{2}, -\frac{m \omega x^{2}}{{{\hbar}}}\right) + K_{2} e^{\left(\frac{m \omega x^{2}}{2 \, {{\hbar}}}\right)} {\rm kummer}_{u}\left(\frac{{{\hbar}} \omega + 2 \, {{\mathcal{E}}}}{4 \, {{\hbar}} \omega}, \frac{1}{2}, -\frac{m \omega x^{2}}{{{\hbar}}}\right)$

The Kummer functions are confluent hypergeometric functions. At present SageMath cannot evaluate these functions so I'll enter the first few wave functions for you manually below.

#Harmonic Oscillator Wavefunction and energy
m,omega=var('m,omega')
hbar=var('hbar',latex_name=r'\hbar')                                  #Declares hbar as a variable and displays hbar symbol when using show() and pretty_print
alpha=m*omega/hbar                               #If alpha is defined _after_ the wave functions are defined, you would need to use .substitute(alpha=m*omega/hbar)

psi_HO0(x) = (alpha/pi)^(1/4)*e^(-alpha*x^2/2)                        #n=0 state of 1D harmonic oscillator
psi_HO1(x) = (4*alpha^3/pi)^(1/4)*x*e^(-alpha*x^2/2)                  #n=1 state of 1D harmonic oscillator
psi_HO2(x) = (alpha/(4*pi))^(1/4)*(2*alpha*x^2-1)*e^(-alpha*x^2/2)    #n=2 state of 1D harmonic oscillator
psi_HO3(x) = (alpha^3/(9*pi))^(1/4)*(2*alpha*x^3-3*x)*e^(-alpha*x^2/2)#n=3 state of 1D harmonic oscillator
psi_HO4(x) = (alpha/(36*pi))^(1/4)*(2*alpha^2*x^4-6*alpha*x^2+3/2)*e^(-alpha*x^2/2)  #n=4 state of 1D harmonic oscillator
psi_HO5(x) = (alpha/(225*pi))^(1/4)*(2*alpha^(5/2)*x^5-10*alpha^(3/2)*x^3+15*alpha^(1/2)*x/2)*e^(-alpha*x^2/2)  #n=5 state of 1D harmonic oscillator


def E_n(n,omega):                                                     #1D harmonic oscillator energy function
    return (n+1/2)*hbar*omega

def H(psi):                                                           #1D harmonic oscillator Hamiltonian
    return -hbar^2/(2*m)*diff(psi,x,2)+1/2*m*omega^2*x^2*psi

In general, the quantum harmonic oscillator wave function can be found from the following relationships:

[\psi_n(x) = \left( \frac{1}{2^n n!} \right) \left( \frac{\alpha}{\pi} \right) H_n(\alpha^{1/2}x) e^{-\alpha x^2/2} ]

where the Hermite polynomial $H_n(x)$ is found from

[H_n(x) = (-1)^n e^{x^2} \frac{d^n}{dx^n} e^{-x^2} ]

For You To Try:

Prove that the quantum harmonic oscillator wave functions are normalized

#Insert your work here

Fancier 2D Plotting

• You can add titles, label axes, change label size, and move the title of plots

• You can specify different line styles:

  • Solid lines '-'

  • Dashed lines '--'

  • Dash dot '-.'

  • Dotted ':'

• If you assign a plot to a variable name, you can add on additional plots and features.

  • In Python, P+=text... is the same as P=P+text.... This is called in-place addition

  • Text can be added to a plot using text('something profound', (x,y), color='red) where x and y are the location of the midpoint of the text.

  • Arrows can be added to point out features by specifying the location of the head and tail of the arrow arrow2d((x1,y1),(x2,y2))

  • Multiple plots can be combined into a single graph by adding plot() commands toP

  Plot the ground state wavefunction of the one-dimensional harmonic oscillator

psi0_plot(x)=psi_HO0(x).substitute(m=1,omega=1,hbar=1)      #To plot a function, all variables must be defined - chose variables = 1 for simplicity

P=plot(psi0_plot(x),[x,-2,2],linestyle='--',title='n=1 State of QHO', axes_labels=['Position x', '$\psi(x)$'],axes_labels_size=1, title_pos=(.3,1))
P+= text('Wave function peak', (1, 1.05), color='red')
P+=arrow2d((1, 1), (0, .75))
P.show()

For You To Try:

• Plot the $n=1$, $n=2$, and $n=3$ wave functions for the quantum harmonic oscillator on the same plot

• Describe how the number of nodes in the wave function relates to the value of n

• Draw arrows to indicate the location of the nodes of the wave function

#Insert your code here

Calculate harmonic oscillator energy from Schrodinger equation and wave function

We can check that the wave functions satisfy Schrodinger's equation and give the correct energy. The equation is

[\hat{H}\psi(x) = E \psi(x).]

Since $E$ is a number (unlike $\hat{H}$, which is an operator), we can divide both sides of the equation by $\psi(x)$ and cancel the terms on the right hand side to get an equation for the energy.

[\frac{\hat{H}\psi(x)}{\psi(x)} = E.]

f=H(psi_HO0(x))/psi_HO0(x)
show(f)

$\displaystyle \frac{{\left(m \omega^{2} x^{2} \left(\frac{m \omega}{\pi {{\hbar}}}\right)^{\frac{1}{4}} e^{\left(-\frac{m \omega x^{2}}{2 \, {{\hbar}}}\right)} - \frac{{\left(\frac{m^{2} \omega^{2} x^{2} \left(\frac{m \omega}{\pi {{\hbar}}}\right)^{\frac{1}{4}} e^{\left(-\frac{m \omega x^{2}}{2 \, {{\hbar}}}\right)}}{{{\hbar}}^{2}} - \frac{m \omega \left(\frac{m \omega}{\pi {{\hbar}}}\right)^{\frac{1}{4}} e^{\left(-\frac{m \omega x^{2}}{2 \, {{\hbar}}}\right)}}{{{\hbar}}}\right)} {{\hbar}}^{2}}{m}\right)} e^{\left(\frac{m \omega x^{2}}{2 \, {{\hbar}}}\right)}}{2 \, \left(\frac{m \omega}{\pi {{\hbar}}}\right)^{\frac{1}{4}}}$

%md

# <img src="sagemathcloud_logo1.png" width="48"> Simplify Symbolic Expressions

- The `.expand()` command will distribute terms and expand powers in an expression
    - Useful for getting terms to cancel out
- The `.canonicalize_radical()` command returns the _canonical_ form of logs and exponents
- Other simplification commands include:
    - `.combine()`,
    - `.simplify_fill()`, 
    - `.simplify_trig()`, 
    - `.reduce_trig()`, 
    - `.simplify_log()`, 
    - `.trig_expand()`

This result does not look like what we would expect.  However, computer algebra systems like Maxima (the heart of SageMath computer algebra) sometimes have issues with non-integer powers.  One trick is to raise this function to the nth power (in this case it would be to the 4th power) to get rid of non-integer powers, simplify the equation, then take the nth root to get the final result.  Let's try this:

Simplify Symbolic Expressions

• The .expand() command will distribute terms and expand powers in an expression

  • Useful for getting terms to cancel out

• The .canonicalize_radical() command returns the canonical form of logs and exponents

• Other simplification commands include:

  • .combine(),

  • .simplify_fill(),

  • .simplify_trig(),

  • .reduce_trig(),

  • .simplify_log(),

  • .trig_expand()

This result does not look like what we would expect. However, computer algebra systems like Maxima (the heart of SageMath computer algebra) sometimes have issues with non-integer powers. One trick is to raise this function to the nth power (in this case it would be to the 4th power) to get rid of non-integer powers, simplify the equation, then take the nth root to get the final result. Let's try this:

f4=(f^4).expand()    #Fourth power

pretty_print('Taking the fourth power of our original result simplifies things down to $f^4=$', f4) 

f1=(f4^(1/4)).canonicalize_radical()    #.canonicalize_radical() simplifies terms involving exponents

pretty_print('Taking the fourth root and using the .canonicalize_radical() command gives us ', f1)
pretty_print('The energy of the n=1 state from the energy function E_n(n,omega) is $E_n=$', E_n(0,omega))

Taking the fourth power of our original result simplifies things down to $f^4=$ $\displaystyle \frac{1}{16} \, {{\hbar}}^{4} \omega^{4}$

Taking the fourth root and using the .canonicalize_radical() command gives us $\displaystyle \frac{1}{2} \, {{\hbar}} \omega$

The energy of the n=1 state from the energy function E_n(n,omega) is $E_n=$ $\displaystyle \frac{1}{2} \, {{\hbar}} \omega$

For you to try:

• Use the Hamiltonian to find the energies of the (n=1), (n=2), and (n=3) states.

  • Remember that (E = \frac{H \psi(x)}{\psi(x)})

  • If your resulting equation is a huge mess, try squaring, expanding, then simplifying the equation

• Check that your results are correct by using the function E_n(n,omega)

#Insert your work here

Calculate the commutator $[\hat{x},\hat{p}]$

We can define the momentum operator $\hat{p}$ using def p:. The operator $\hat{x}$ is just multiplication by the variable $x$ so we don't need to define a separate operator. One thing to keep in mind is that we must explicitly apply the operators (i.e. All terms that the operator operates on must be inside the parentheses as an argument to the operator function)

def p(psi):   #momentum operator
    return -i*hbar*diff(psi,x)


psi=function('psi')(x)   #Dummy function


com=x*p(psi)-p(x*psi)   #Calculate the commutation relation [x,p]
pretty_print('The commutation relation between $\hat{x}$ and $\hat{p}$ is $[\hat{x},\hat{p}]\psi(x) = $', com.expand())  #.expand() required to simplify expression

The commutation relation between $\hat{x}$ and $\hat{p}$ is $[\hat{x},\hat{p}]\psi(x) =$ $\displaystyle i \, {{\hbar}} \psi\left(x\right)$

Raising and Lowering Operators

The raising and lowering operators are powerful tools for exploring the quantum harmonic oscillator. Most quantum physics texts will cover these in detail. The following cells will explore a few things you can do in SageMath with the raising and lowering operators. The cell below defines the operators.

Note: The output from the cell below includes a term $D[0]\psi(x)$. This is SageMath notation for the derivative with respect to the first variable in the function (remember that Python starts numbering at zero). If $\psi$ was a function of $x$ and $y$, then the derivative of $\psi(x,y)$ with respect to $y$ would be denoted $D[1]\psi(x,y)$.

psi,m,L,omega, f=var('psi,m,L,omega,f')                          #Declare wave function, mass, and box length as variables

Ener=var('Ener',latex_name=r'\mathcal{E}')                       #Declare Ener as energy variable but will display a script E when using show() and pretty_print
hbar=var('hbar',latex_name=r'\hbar')                             #Declares hbar as a variable and displays hbar symbol when using show() and pretty_print

psi=function('psi')(x)    #Define dummy function

def p(psi):   #momentum operator
    return -i*hbar*diff(psi,x)


def a_plus(psi):   #raising operator
    return 1/sqrt(2*hbar*m*omega)*(-i*p(psi)+m*omega*x*psi)


def a_minus(psi):   #lowering operator
    return 1/sqrt(2*hbar*m*omega)*(+i*p(psi)+m*omega*x*psi)

show(a_plus(psi))

$\displaystyle \frac{\sqrt{2} {\left(m \omega x \psi\left(x\right) - {{\hbar}} D[0]\left(\psi\right)\left(x\right)\right)}}{2 \, \sqrt{{{\hbar}} m \omega}}$

Demonstrate that $a_+ \psi_{0}(x)= \psi_{1}(x)$

• Any easy way to show two equations are equal is to divide one equation by the other and show the result is equal to 1

  • This cuts down on manipulation of equations to put them in the same format

• To get all terms to cancel we will raise everything to the 4th power and then take the 4th root

• The .canonicalize_radical() command is used to expand the exponent terms

Show that $\frac{a_+ \psi_0(x)}{\psi_1(x)} = 1$

f=(a_plus(psi_HO0(x))/psi_HO1(x))

pretty_print(' \( a_+ \psi_0(x)/\psi_1(x) = \)' ,f)   #It's not clear this is equal to one but ...
f4=(f^4)                                              #... with a little nudge...
f1=(f4^(1/4)).canonicalize_radical()
pretty_print(' \( a_+ \psi_0(x)/\psi_1(x) = \)' ,f1)  #...everything works out nicely

$a_+ \psi_0(x)/\psi_1(x) =$ $\displaystyle \frac{4^{\frac{3}{4}} \sqrt{2} m \omega \left(\frac{m \omega}{\pi {{\hbar}}}\right)^{\frac{1}{4}}}{4 \, \sqrt{{{\hbar}} m \omega} \left(\frac{m^{3} \omega^{3}}{\pi {{\hbar}}^{3}}\right)^{\frac{1}{4}}}$

$a_+ \psi_0(x)/\psi_1(x) =$ $\displaystyle 1$

Demonstrate that $a_+ \psi_{1}(x)= \sqrt{2}\psi_{2}(x)$

• In addition to .canonicalize_radical() we also have to use .combine() to collect terms over a common denominator, and .expand().

Show that $\frac{a_+ \psi_1(x)}{\psi_2(x)} = \sqrt{2}$

f=(a_plus(psi_HO1(x))/psi_HO2(x)).expand().combine()

pretty_print(' \( a_+ \psi_1(x)/\psi_2(x) = \)' ,f) 

f4=(f^4).expand()
f1=(f4^(1/4)).canonicalize_radical()

pretty_print(' \( a_+ \psi_1(x)/\psi_2(x) = \)' ,f1)
pretty_print(' or equivalently \( a_+ \psi_1(x) =\)',f1,' \(\psi_2(x)\)')

$a_+ \psi_1(x)/\psi_2(x) =$ $\displaystyle \frac{2 \, \sqrt{2} \left(\frac{m^{3} \omega^{3}}{{{\hbar}}^{3}}\right)^{\frac{1}{4}} m \omega x^{2} - \sqrt{2} \left(\frac{m^{3} \omega^{3}}{{{\hbar}}^{3}}\right)^{\frac{1}{4}} {{\hbar}}}{\sqrt{{{\hbar}} m \omega} {\left(\frac{2 \, m \omega x^{2}}{{{\hbar}}} - 1\right)} \left(\frac{m \omega}{{{\hbar}}}\right)^{\frac{1}{4}}}$

$a_+ \psi_1(x)/\psi_2(x) =$ $\displaystyle \sqrt{2}$

or equivalently $a_+ \psi_1(x) =$ $\displaystyle \sqrt{2}$ $\psi_2(x)$

For You To Try

Show that Show that $\frac{a_+ \psi_2(x)}{\psi_1(x)} = \sqrt{3}$

#Insert your work here

Projects and Problem Ideas

• Use $a_+$ to generate higher quantum harmonic oscillator wave functions and then use $H=\hbar \omega(a_- a_+ - \frac{1}{2})$ to determine the energy of those states.

• Rewrite $\hat{x}$ and $\hat{p}$ in terms of $a_-$ and $a_+$ and use the raising and lowering operators to calculate the expectation values $<x>$ and $<p>$. Compare your results to what you get from $\int \psi^*(x) x \psi(x) dx$ and $\int \psi^*(x) p \psi(x) dx$.