iTranslated by AI

The content below is an AI-generated translation. This is an experimental feature, and may contain errors. View original article

🎸

Visualization with Julia: Anharmonic Oscillator

2022/01/07に公開

2023/12/18

Visualize the solutions to the Schrödinger equation for a particle moving in a Morse potential.

Packages

You need to install the packages beforehand by running add Plots, add LaTeXStrings, add SpecialPolynomials, and add SpecialFunctions in Julia's package mode. Also, declare using Plots, using LaTeXStrings, using SpecialPolynomials, and using SpecialFunctions in your notebook. For instructions on installing and using Jupyter Notebook, see this article, and for how to use Plots.jl, see this article.

Packages

# using Pkg
# Pkg.add("Plots")
# Pkg.add("LaTeXStrings")
# Pkg.add("SpecialPolynomials")
# Pkg.add("SpecialFunctions")
using Plots
using LaTeXStrings
using SpecialPolynomials
using SpecialFunctions

Time-Independent Schrödinger Equation

Morse Potential

The Morse potential, with symbols slightly modified from Equation (4) in Morse (1929), is given by

V(r) = D_\mathrm{e} \left( \mathrm{e}^{-2a(r-r_e)} - 2\mathrm{e}^{-a(r-r_e)} \right)

and the Hamiltonian for a particle of mass $\mu$ moving in this potential is

\hat{H} = - \frac{\hbar^2}{2\mu} \frac{\mathrm{d}^2}{\mathrm{d}r ^2} + V(r)

The solution to the corresponding time-independent Schrödinger equation

\hat{H}\psi(r) = E \psi(r)

is given by the following, modifying Toda (1999) and Wikipedia respectively:

\psi_n(\xi) = N_n z^{\lambda-n-1/2} \mathrm{e}^{-z/2} L_n^{(2\lambda-2n-1)}(\xi)\\ E_n = -D_\mathrm{e} + \hbar \omega \left( n + \frac{1}{2} \right) - \chi \hbar \omega \left( n + \frac{1}{2} \right)^2

where $\xi=2\lambda\mathrm{e}^{-a(r-r_e)}$ , $\omega=\sqrt{k/µ}$ , $k=2D_\mathrm{e}a^2$ , $\lambda=\frac{\sqrt{2mD_\mathrm{e}}}{a\hbar}$ , $\chi=\frac{\hbar\omega}{4D_\mathrm{e}}$ , $N_n=\sqrt{\frac{n!(2\lambda-2n-1)a}{\Gamma(2\lambda-n)}}$ , and $L_n^{(\alpha)}$ is the generalized Laguerre polynomial. Note that these are not the associated Laguerre polynomials (see here for details). Additionally, I added $a$ to the normalization constant $N_n$ as it seemed to be missing in Wikipedia.

Potential

V(r, re, De, a) = De*( exp(-2*a*(r-re)) -2*exp(-a*(r-re)) )

Eigenfunctions

# using SpecialFunctions
# using SpecialPolynomials
Γ(z) = try; gamma(z); catch; Inf; end

function ψ(n, r, re ,De, a, µ)
    ℏ  = 1.0
    λ = sqrt(2*µ*De)/(a*ℏ)
    ξ = 2*λ*exp(-a*(r-re))                                    
    s  = 2*λ-2*n-1
    N  = sqrt(factorial(n)s*a/Γ(s+n+1))
    return N*ξ^(s/2)*exp(-ξ/2)*basis(Laguerre{s}, n)(ξ)
end

Eigenvalues

function E(n, De, a, µ)
    ℏ = 1.0
    k = 2*a^2*De
    ω = sqrt(k/µ)
    χ = ℏ*ω/(4*De)
    return - De + ℏ*ω*(n+1/2) - χ*ℏ*ω*(n+1/2)^2
end

Parameters

re = 1.38436
De = 1.13133356259 - 1
k = 12.772/15.569141/2
a = sqrt(k/(2*De))
µ = 1/(1/1837.108542 + 1/1837.108542)

println("re = ", re)
println("De = ", De)
println("k  = ", k)
println("a  = ", a)
println("µ  = ", µ)

Output

re = 1.38436
De = 0.13133356259000006
k  = 0.4101703491541377
a  = 1.249623750214457
µ  = 918.554271

Using the functions and parameters above, the Morse potential can be depicted as follows. It can be seen that it asymptotically approaches $V(\infty)=0$ .

# Axis labels, plot range, legend position
plot(xlabel=L"r", ylabel=L"V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# Potential
plot!(-3:0.02:7, r->V(r,re,De,a), label="", lc="#578FC7", lw=2)

# Display
plot!() |> display

In the Morse potential, the condition $E_{n+1} - E_{n} > 0$ fails for levels that are not finite bound states. Listing these shows that there are bound states from $n=0$ to $11$ .

for n in 0:15
    println(n, "\t", E(n+1,De,a,µ)-E(n,De,a,µ) > 0)
end

0	true
1	true
2	true
3	true
4	true
5	true
6	true
7	true
8	true
9	true
10	true
11	true
12	false
13	false
14	false
15	false

Depicting all the energy levels of the bound states:

# Axis labels, plot range, legend position
plot(xlabel=L"r", ylabel=L"V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# Potential
plot!(0:0.02:7, r->V(r,re,De,a), label="", lc="#578FC7", lw=2)

# Energy
for n in 11:-1:0
    plot!(0:0.02:7, r-> E(n,De,a,µ)<V(r,re,De,a) ? NaN : E(n,De,a,µ), label="", lc="#578FC7", lw=1)
end

# Display
plot!() |> display

Shift the wavefunctions to the height of their energy levels, scale them vertically, and add them to the plot. Here, $0.004\psi(x)+E_n$ is plotted. It can be seen that the quantum number $n$ corresponds to the number of nodes.

# Axis labels, plot range, legend position
plot(xlabel=L"r", ylabel=L"V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# Potential
plot!(-3:0.02:7, r->V(r,re,De,a), label="", lc="#000000", lw=2)

for n in 11:-1:0
    # Energy (dashed line)
    hline!([E(n,De,a,µ)], label="", lc="#000000", ls=:dash)
    # Energy (solid line)
    plot!(0:0.02:7, r-> E(n,De,a,µ)<V(r,re,De,a) ? NaN : E(n,De,a,µ), label="", lc="#000000", lw=2)
    # Wavefunction
    plot!(0:0.02:7, r->0.004*ψ(n,r,re,De,a,µ)+E(n,De,a,µ), label=latexstring("n=$n"), lc=n+1, lw=2)
end

# Display
plot!() |> display

Verification by Numerical Integration

Since it is often impossible to judge whether the eigenfunctions are correctly declared just by looking at them, we perform a simple calculation to check if the normalization condition $\int |\psi_n(x)|^2 \mathrm{d}x=1$ is satisfied. This also helps to verify if the library is used correctly, if there are any bugs, or if there are errors in the literature. The integral is discretized as follows. The range of $i$ is taken to be sufficiently wide. Since this is numerical integration, it won't be exactly 1, but if it is roughly 1, it should be fine.

\int_0^\infty |\psi(r)|^2 \mathrm{d}r \simeq \sum_{i=0}^\infty |\psi(i \Delta r)|^2 \Delta r

Numerical integration routine

# Numerical integration routine
integral(func, x_min, x_max, dx) = sum(func.(x_min:dx:x_max)*dx)

# Verification of normalization condition

println("n \t∫|ψ(r)|²dr")
for n in 0:11
    println(n, "\t", integral(r-\u003eψ(n,r,re,De,a,µ)^2, -10, 50, 0.001))
end

n 	∫|ψ(r)|²dr
0	1.0000000000000002
1	1.0
2	0.9999999999999998
3	1.0000000000000004
4	1.0
5	0.9999999999999997
6	0.9999999999999999
7	0.9999999999999997
8	0.9999999999999989
9	1.0
10	0.9999999999999999
11	0.9999999999999988

Next, we evaluate the energy through numerical integration. Unlike the harmonic oscillator, the virial theorem cannot be applied here, so we must deal with the kinetic term directly. The second derivative can be approximated using the finite difference method as:

\begin{align} \frac{\mathrm{d}^2\psi(r)}{\mathrm{d}r ^2} \u0026= \frac{ \frac{\psi(r+\Delta r) - \psi(r)}{\Delta r} - \frac{\psi(r) - \psi(r-\Delta r)}{\Delta r} }{\Delta r} \\ \u0026= \frac{\psi(r+\Delta r) - 2\psi(r) + \psi(r-\Delta r)}{\Delta r ^2} \end{align}

Using the fact that the wavefunction is a real function, $\psi^*(r)=\psi(r)$ , we can discretize the expression as:

\begin{align} \int_0^\infty \psi^*(r) \hat{H} \psi(r) \mathrm{d}r \u0026= \int_0^\infty \psi^*(r) \left[ - \frac{\hbar^2}{2\mu} \frac{\mathrm{d}^2}{\mathrm{d}r ^2} + V(r) \right] \psi(r) \mathrm{d}r \\ \u0026= - \frac{\hbar^2}{2\mu} \int_0^\infty \psi(r) \frac{\mathrm{d}^2 \psi(r)}{\mathrm{d}r ^2} \mathrm{d}r + \int_0^\infty V(r) |\psi(r)|^2 \mathrm{d}r \\ \u0026= - \frac{\hbar^2}{2\mu} \sum_{i=1}^\infty \psi(i\Delta r)\frac{\psi(i\Delta r+\Delta r) - 2\psi(i\Delta r) + \psi(i\Delta r-\Delta r)}{\Delta r ^2} \Delta r + \sum_{i=1}^\infty V(i\Delta r) |\psi(i\Delta r)|^2 \Delta r \end{align}

The starting index was changed to $i=1$ to ensure that the argument of $\psi(r)$ does not become negative. Since $\Delta r \ll 1$ , there should be little difference between $i=0$ and $i=1$ . Depending on the precision of $\Delta r$ , it can be verified that several digits match.

println("n \tE\t\t\t∫ψHψdr")

dr = 0.005
ψVψ(n,r,re,De,a,µ)    = V(r,re,De,a) * ψ(n,r,re,De,a,µ)^2
ψTψ(n,r,re,De,a,µ,dr) = -1/(2*µ) * ψ(n,r,re,De,a,µ) * (ψ(n,r+dr,re,De,a,µ)-2*ψ(n,r,re,De,a,µ)+ψ(n,r-dr,re,De,a,µ))/dr^2
ψHψ(n,r,re,De,a,µ,dr) = ψTψ(n,r,re,De,a,µ,dr) + ψVψ(n,r,re,De,a,µ)

for n in 0:11
    println(n, "\t", E(n,De,a,µ), "\t", integral(r-\u003eψHψ(n,r,re,De,a,µ,dr), dr, 100, dr))
end

n 	E			∫ψHψdr
0	-0.12098032966578513	-0.12098063308527979
1	-0.10154887790721412	-0.10155018470683336
2	-0.08381744493512158	-0.08382034362426742
3	-0.06778603074950751	-0.06779072698018129
4	-0.05345463535037193	-0.053461018312416515
5	-0.040823258737714835	-0.040830967539695506
6	-0.0298919009115362	-0.02990039094983932
7	-0.020660561871836063	-0.020669171190581332
8	-0.013129241618614396	-0.0131372572629767
9	-0.007297940151871213	-0.0073046645173990635
10	-0.003166657471606507	-0.0031714746521237978
11	-0.0007353935778202786	-0.0007378357145039691

Harmonic Approximation

The harmonic oscillator potential corresponding to the harmonic approximation of this potential, along with the eigenvalues and eigenfunctions of its Hamiltonian, are given by

V(r) = V_0 - D_\mathrm{e} + \frac{1}{2} k (r - r_\mathrm{e})^2\\ \psi_n(r) = A_n H_n(\xi) \exp{\left( -\frac{\xi^2}{2} \right)} \\ E_n = V_0 -D_\mathrm{e} + \hbar \omega \left( n + \frac{1}{2} \right)

where $\omega:=\sqrt{k/m}$ , $\xi=\sqrt{\frac{m\omega}{\hbar}} (r - r_\mathrm{e})$ , $A_n=\sqrt{\frac{1}{n! 2^n} \sqrt{\frac{m\omega}{\pi\hbar}}}$ , and $H_n$ is the Hermite polynomial. Also,

a = \sqrt{k/2D_\mathrm{e}}

Comparing the energy with the Morse potential, it can be seen that the difference lies only in the third term (the anharmonic term). Additionally, although it is called the harmonic approximation, it is often more practical to determine the spring constant $k$ first and then determine $a$ using the above equation.

module HO
using SpecialPolynomials

# Potential
V(r, re, De, k) = - De + 1/2*k*(r-re)^2

# Eigenfunctions
# using SpecialPolynomials
function ψ(r,re,n,k,m)
    ℏ = 1.0
    ω = sqrt(k/m) # mu = 1 / (1/m1 + 1/m2)
    A = sqrt(1/(factorial(n)*2^n)*sqrt(m*ω/(π*ℏ)))
    ξ = sqrt(m*ω/ℏ) * (r-re)
    return A*basis(Hermite, n)(ξ)*exp(-ξ^2/2)
end

# Eigenvalues
function E(n,De,k,m)
    ℏ = 1
    ω = sqrt(k/m)
    return - De + ℏ*ω*(n+1/2)
end

end

Depicting it and comparing it with the Morse potential, we can see that they have a very similar shape near the minimum point.

# Axis labels, plot range, legend position
plot(xlabel=L"r", ylabel=L"V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# Potential
plot!(0:0.02:7, r-\u003eV(r,re,De,a), label="Morse Potential", lc="#578FC7", lw=2)
plot!(0:0.02:7, r-\u003eHO.V(r,re,De,k), label="Harmonic Oscillator", lc="#BC1C5F", lw=2)

# Display
plot!() |\u003e display

Next, we plot the energy levels. At lower levels, they match well compared to higher levels. Since the level spacing is constant, the number of bound states is smaller.

# Axis labels, plot range, legend position
plot(xlabel=L"r", ylabel=L"V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# Potential
plot!(0:0.02:7, r-\u003eV(r,re,De,a), label="Morse Potential", lc="#578FC7", lw=2)
plot!(0:0.02:7, r-\u003eHO.V(r,re,De,k), label="Harmonic Oscillator", lc="#BC1C5F", lw=2)

# Energy
for n in 11:-1:0
    plot!(0:0.02:7, r-\u003e E(n,De,a,µ)\u003cV(r,re,De,a) ? NaN : E(n,De,a,µ), label="", lc="#578FC7")
end

for n in 5:-1:0
    plot!(0:0.02:7, r-\u003e HO.E(n,De,k,µ)\u003cHO.V(r,re,De,k) ? NaN : HO.E(n,De,k,µ), label="", lc="#BC1C5F")
end

# Display
plot!() |\u003e display

Next, we plot the wavefunctions of the harmonic oscillator.

# Axis labels, plot range, legend position
plot(xlabel=L"r", ylabel=L"V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# Potential
plot!(-3:0.02:17, r-\u003eV(r,re,De,a), label="", lc="#737B85", lw=1, ls=:dash)
plot!(-3:0.02:17, r-\u003eHO.V(r,re,De,k), label="", lc="#000000", lw=2)

for n in 11:-1:0
    # Energy (solid line)
    plot!(0:0.02:7, r-\u003e E(n,De,a,µ)\u003cV(r,re,De,a) ? NaN : E(n,De,a,µ), label="", lc="#737B85", lw=1, ls=:dash)
end

for n in 5:-1:0
    # Energy (dashed line)
    # hline!([HO.E(n,De,k,µ)], label="", lc="#000000", ls=:dash)
    # Energy (solid line)
    plot!(0:0.02:7, r-\u003e HO.E(n,De,k,µ)\u003cHO.V(r,re,De,k) ? NaN : HO.E(n,De,k,µ), label="", lc="#000000", lw=2)
    # Wavefunction
    plot!(0:0.02:7, r-\u003e0.004*HO.ψ(r,re,n,k,µ)+HO.E(n,De,k,µ), label=latexstring("n=$n"), lc=n+1, lw=2)
end

# Display
plot!() |\u003e display

Comparing the wavefunctions of the harmonic oscillator (dashed lines) and the Morse oscillator (solid lines) yields the following. At lower levels, the anharmonicity is small, and the description by the harmonic approximation is not bad compared to higher levels. At higher vibrational levels, in addition to the energy difference, the tendency for the maximum value of the wavefunction to shift toward the long-range side is also clearly apparent.

# Axis labels, plot range, legend position
plot(xlabel=L"r", ylabel=L"V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# Potential
plot!(-3:0.02:17, r-\u003eV(r,re,De,a), label="", lc="#000000", lw=2)
plot!(-3:0.02:17, r-\u003eHO.V(r,re,De,k), label="", lc="#000000", lw=1, ls=:dot)

# Wavefunction
for n in 5:-1:0
    plot!(0:0.02:7, r-\u003e0.004*HO.ψ(r,re,n,k,µ)+HO.E(n,De,k,µ), label="", lc=n+1, lw=1, ls=:dot)
    plot!(0:0.02:7, r-\u003e0.004*ψ(n,r,re,De,a,µ)+E(n,De,a,µ), label=latexstring("n=$n"), lc=n+1, lw=2)
end

# Display
plot!() |\u003e display

Time-Dependent Schrödinger Equation

The solution to the time-dependent Schrödinger equation

i\hbar \frac{\mathrm{d}}{\mathrm{d}t} \varPsi(x,t) = \hat{H} \varPsi(x,t)

is generally given by

\varPsi(x,t) = \sum_{n=0}^\infty c_n \psi_n(x) \exp\left(-i \frac{E_n}{\hbar} t \right)

That this is a solution can be shown simply by substituting it (verify the conditions under which the summation and the differential operator are interchangeable). Here, $c_n$ are the coefficients defining the initial condition

\varPsi(x,0) = \sum_{n=0}^\infty c_n \psi_n(x)

Time-dependent wavefunction

function Ψ(N,t,r,re,De,a,µ) # ← Note that this is an uppercase Ψ
    sum = 0.0 + im*0.0
    for n in 0:N
        sum += c(n) * ψ(n,r,re,De,a,µ) * exp(im*E(n,De,a,µ)*t)
    end
    return sum
end

# Depicting time evolution
function draw(num)
    anim = Animation()
    N = 10
    for i in 0:100
        # Time
        t = 10.0*i
        # Plot range, title, etc.
        plt = plot(title="t=$t", xlabel="r", ylabel="V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))
        # Potential
        plot!(0:0.01:17, r->V(r,re,De,a), label="", lc="#000000", lw=2)
        # Real part Re Ψ(x,t)
        plot!(plt, 0:0.01:7, r->real(Ψ(N,t,r,re,De,a,µ))*0.05 + sum(c.(0:N) .* E.(0:N,De,a,µ)),  label="Re", lc=1)
        # Imaginary part Im Ψ(x,t)
        plot!(plt, 0:0.01:7, r->imag(Ψ(N,t,r,re,De,a,µ))*0.05 + sum(c.(0:N) .* E.(0:N,De,a,µ)),  label="Im", lc=2)
        # Probability density |Ψ(x,t)|^2
        plot!(plt, 0:0.01:7, r->abs(Ψ(N,t,r,re,De,a,µ))^2*0.05 + sum(c.(0:N) .* E.(0:N,De,a,µ)), label="P(x,t)", lc=3, lw=2)
        # Add frame
        frame(anim, plt)
    end
    gif(anim, "TDSE_$(num).gif", fps = 20)
end

First, if the initial condition is purely the ground state ( $c_0=1, c_n=0 (n>0)$ ), the real and imaginary parts alternate, and the probability density does not change over time. This can be easily understood by studying the finite difference method.

c(n) = n==0 ? 1 : 0
draw(1)

Next, considering an initial condition where the ground state and the first excited state are mixed in equal proportions ( $c_0=c_1=0.5, c_n=0 (n>1)$ ), it begins to exhibit oscillatory behavior.

c(n) = n<2 ? 0.5 : 0
draw(2)

Finally, if we set $c_0=\frac{1}{2}, c_1=\frac{1}{4}, c_2=\frac{1}{8}, \cdots, c_n=\frac{1}{2^{n+1}}$ , we can observe a strange motion that differs from that of a harmonic oscillator.

c(n) = 2.0^(-n-1)
draw(3)

Operating Environment

versioninfo()

Output

    Julia Version 1.6.2
    Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
    Platform Info:
      OS: Windows (x86_64-w64-mingw32)
      CPU: Intel(R) Core(TM) i7-4650U CPU @ 1.70GHz
      WORD_SIZE: 64
      LIBM: libopenlibm
      LLVM: libLLVM-11.0.1 (ORCJIT, haswell)

References

The Jupyter Notebook data that this article is based on can be found at the following link.

iTranslated by AI

Visualization with Julia: Anharmonic Oscillator

Packages

Time-Independent Schrödinger Equation

Morse Potential

Verification by Numerical Integration

Harmonic Approximation

Time-Dependent Schrödinger Equation

Operating Environment

References

Julia Visualization Series

Discussion