Correct Setup for Finite Difference to Calculate Quantum Tunneling

[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 probability 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?

 

[DrClaude]

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.)

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.

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.

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.

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.

However, when I recompute for different energies, the tunneling probability 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.

 

[Tertius]

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.