# Correct Setup for Finite Difference to Calculate Quantum Tunneling

• Tertius
Tertius
Homework Statement
I am discretizing the Schrodinger equation and solving for a quantum tunneling probability
Relevant Equations
Schrodinger equation
I thought I solved the problem in answering my own post a few days ago, but the tunneling probability vs. energy trend is clearly wrong. I've remade the post because I have totally changed my approach and need a better understanding of the boundary setup.

Overall description: a plane wave approaches a barrier, tunnels, and a plane wave exits. I am doing a full 1D simulation so I can eventually use non-standard potentials, but I'm starting with a square potential.

Setup of Solution:
$$-\frac{\hbar^{2}}{2m} \frac{d^{2}}{dx^{2}} \psi(x) +V(x)\psi(x)=E\psi(x)$$
$$H=-\frac{\hbar^{2}}{2m}+V(x)$$
The tri-diagonal matrix setup after applying finite derivates is ##\frac{h^2}{2ma^2}+V_i## for diagonals and ##\frac{-h^2}{2ma^2}## for the off-diagonals.

I have not included any boundary conditions explicitly. The standard setup for this problem is $$H\psi=E\psi$$ where H is the Hamiltonian and E are the energy eigenvalues. I don't need to compute the eigenvalues, however because this is an unbounded system. The energy eigenvalue for this system should be the kinetic energy of the incoming particle. This allows me to set it up as $$(H-IE)\psi=b$$ The vector b is used for the right side as an empty vector. ##I## is the identity matrix.

When the vector ##b## is actually all 0, the solution is always trivial. I have to include a very small, non-zero value somewhere in this vector for it to produce a solution that looks like the attached picture. However, when I recompute for different energies, the tunneling probablity changes in a weird, oscillatory way as shown in the other attached graph.

Question:
What is wrong with the assumptions of this setup? Do I need to include boundaries specifically, even though the solutions are plane waves?

Last edited:
Tertius said:
Overall description: a plane wave approaches a barrier, tunnels, and a plane wave exits. I am doing a full 1D simulation so I can eventually use non-standard potentials, but I'm starting with a square potential.
A plane wave cannot "approach" a barrier. A plane wave is a solution that covers all space. Do you mean a wave packet?

(Note: I know that when dealing with scattering, the incoming state is taken as a plane wave, but that is in the limit ##t \rightarrow \infty##, which is not suitable for a numerical simulation.)

Tertius said:
Setup of Solution:
$$-\frac{\hbar^{2}}{2m} \frac{d^{2}}{dx^{2}} \psi(x) +V(x)\psi(x)=E\psi(x)$$
$$H=-\frac{\hbar^{2}}{2m}+V(x)$$
The term ##\frac{d^{2}}{dx^{2}}## is missing here.

Tertius said:
The tri-diagonal matrix setup after applying finite derivates is ##\frac{h^2}{2ma^2}+V_i## for diagonals and ##\frac{-h^2}{2ma^2}## for the off-diagonals.
The kinetic term on the diagonal is off by a factor of 2.

Tertius said:
The standard setup for this problem is $$H\psi=E\psi$$ where H is the Hamiltonian and E are the energy eigenvalues.
Depends what you are after. The time-dependent Schrödinger equation might be more appropriate.

Tertius said:
I don't need to compute the eigenvalues, however because this is an unbounded system. The energy eigenvalue for this system should be the kinetic energy of the incoming particle.
If you solve the time-independent Schrödinger equation as above, with the potential as part of the Hamiltonian, you will get energy eigenstates that are not free-particle eigenstates, hence the energy won't correspond to the kinetic energy.

Tertius said:
However, when I recompute for different energies, the tunneling probablity changes in a weird, oscillatory way as shown in the other attached graph.
How do you calculate that probability?

Question:
What is wrong with the assumptions of this setup? Do I need to include boundaries specifically, even though the solutions are plane waves?
[/QUOTE]
Again, the solutions are not plane waves. AS set up, the problem assumes hard walls at the extremities of the grid, so the eigenfunctions should go to zero at both ends. You can solve the full eigenvalue problem ##H\psi=E\ps## and get the corresponding discrete set of eigenstates.

WWGD
DrClaude said:
A plane wave cannot "approach" a barrier. A plane wave is a solution that covers all space. Do you mean a wave packet?

(Note: I know that when dealing with scattering, the incoming state is taken as a plane wave, but that is in the limit ##t \rightarrow \infty##, which is not suitable for a numerical simulation.)

The term ##\frac{d^{2}}{dx^{2}}## is missing here.

The kinetic term on the diagonal is off by a factor of 2.

Depends what you are after. The time-dependent Schrödinger equation might be more appropriate.

If you solve the time-independent Schrödinger equation as above, with the potential as part of the Hamiltonian, you will get energy eigenstates that are not free-particle eigenstates, hence the energy won't correspond to the kinetic energy.

How do you calculate that probability?

Question:
What is wrong with the assumptions of this setup? Do I need to include boundaries specifically, even though the solutions are plane waves?
Again, the solutions are not plane waves. AS set up, the problem assumes hard walls at the extremities of the grid, so the eigenfunctions should go to zero at both ends. You can solve the full eigenvalue problem ##H\psi=E\ps## and get the corresponding discrete set of eigenstates.
[/QUOTE]

Thanks for the response.

Ok so trying to use plane waves in a numerical simulation necessitates a boundary, which means they aren't plane waves. So, the current setup was really just solving the equation between two boundaries with a potential in the middle. And it would require solving for the eigensystem and the eigenstates may not match the kinetic energy of the incoming particle.

I am open to changing my approach, my intention is to just have accurate numerical estimates of tunneling probabilities so I can eventually use non-standard potentials. I am reading through https://www.reed.edu/physics/faculty/wheeler/documents/Quantum Mechanics/Miscellaneous Essays/Gaussian Wavepackets.pdf

and I have a Gaussian wave packet simulation built, and it appears more robust.

I suppose my question now is, how to compute the tunneling probability from numerical results. Do I simply compute ##\psi^* \psi## at the front and back of the barrier (and take the ratio)? Basically I would just compute that value at each point in time in the simulation and see where is it maximal.

Last edited:

• Introductory Physics Homework Help
Replies
3
Views
911
• Introductory Physics Homework Help
Replies
2
Views
639
• Introductory Physics Homework Help
Replies
8
Views
1K
• Introductory Physics Homework Help
Replies
9
Views
859
• Introductory Physics Homework Help
Replies
22
Views
2K
Replies
1
Views
338
• Introductory Physics Homework Help
Replies
4
Views
1K
• Introductory Physics Homework Help
Replies
3
Views
1K
• Quantum Physics
Replies
7
Views
748
• Introductory Physics Homework Help
Replies
8
Views
5K