Applied QM - Gaussian Potential

We are looking for solutions of the differential equation $-y''(x) + e^{-x^2}y(x) \; = \; E y(x) \hspace{2cm} 0 \le x \le 5$ with $y(0) = 1$ , $y'(0) = 0$ (the ground state will be an even function) and $y(5) = 0$ .

Solving the differential equation, with the initial conditions $y(0)=1$ , $y'(0)=0$ , numerically, we can solve numerically to find the smallest positive value of $E$ that gives $y(5)=0$ ; this is $E = 0.331874071$ .

For this value of $E$ , we then calculate $\frac{\int_0^1 y(x)^2\,dx}{\int_0^5 y(x)^2\,dx} \; = \; \boxed{0.157447}$ for the required probability.

We have already seen some nice solutions by @Mark Hennings and @Karan Chatrath. The only thing I will add here is some plots for the higher-energy modes. The energies are as follows:

$E_1 \approx 0.332 \\ E_2 \approx 0.445 \\ E_3 \approx 1.189 \\ E_4 \approx 1.709 \\ E_5 \approx 2.684 \\ E_6 \approx 3.719$

The first two modes "prefer" the low-potential regions, analogously with classical mechanics. From then on, the odd-numbered modes have a local maximum in the probability density coinciding with the potential peak, and the even-numbered modes do not.

Post-script :

The energy levels above were found by imposing the boundary conditions $\Psi(-5) = 0$ and $\frac{d}{dx} \Psi (-5) \neq 0$ . If we instead impose the conditions $\Psi(0) = 0$ and $\frac{d}{dx} \Psi (0) \neq 0$ , we get another set of energies, some of which are basically in common with the first set. Notice that $E_1$ , $E_2$ , and $E_3$ from this set correspond to $E_2$ , $E_4$ , and $E_6$ from the other set. The probability density plots for these modes are given below. Notice that all of them now "avoid" the high potential region.

$E_1 \approx 0.445 \\ E_2 \approx 1.712 \\ E_3 \approx 3.725 \\ E_4 \approx 6.495 \\ E_5 \approx 10.048 \\ E_6 \approx 14.389$

I have nothing new to add to the solution posted by @Mark Hennings, except that this is the first time I have used the Brilliant coding environment.

The differential equation can be re-written as a system of first order ODE's by taking $y_1(x) = \psi(x)$ and $y_2(x) = \dot{\psi}(x)$ :

$\dot{y}_1 = y_2$ $\dot{y}_2 = (e^{-x^2} - E)y_1$

The dot notation denotes derivative with respect to $x$ .

$\left[\begin{matrix}\dot{y}_1 \\ \dot{y}_2\end{matrix}\right] = \left[\begin{matrix}y_2 \\ (e^{-x^2} - E)y_1\end{matrix}\right]$ $\dot{y} = f(x,y)$

It has been established in previous problems that the shape of the wave function is independent of the spatial derivative. So taking boundary conditions at $x=-5$ as:

$y(-5) = \left[\begin{matrix}0\\1\end{matrix}\right]$

The system of ODE's is numerically integrated using the Euler step as such:

$y(k+1) = y(k) + f(x(k),y(k))\delta x$

As $x$ varies from -5 to 5. The value of energy $E$ was incrementally increased until $\psi(5) = 0$ . This led to:

$E = 0.33185$

The probability is computed as such:

$P = \frac{\int_{-1}^{1}\psi^2(x)dx}{\int_{-5}^{5}\psi^2(x)dx}$

As for the bonus question, the wave function indicates that the particle is unlikely to be found near the area of maximum PE. This is because the region around the maximum corresponds to an unstable equilibrium. This is in line with classical mechanics.

Simulation code is as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

import math
import numpy as np

#constants
E = 0.33185

# Initialize x

x = -5.0
dx = 10.0**(-5.0)
xs = np.array([[0.0],[1.0]])
A1 = 0
A2 = 0

# Numerical Integration - Euler step:
while x <= 5.0:
    xs = xs +  np.array([xs[1],xs[0]*(math.exp(-x**2)-E)])* dx

    # Computation of area under the wave-function squared:
    A1 = A1 + xs[0]**2*dx

    # Computation of area under the wave function squared when x is between -1 and 1:
    if x >= -1 and x <= 1:
        A2 = A2 + (xs[0]**2)*dx

    x = x + dx

P = A2/A1

print(P)