Best method for numerically solving Schrodinger equation for quantum tunneling w/ arbitrary potential

Question

Goal: to compute accurate tunneling probabilities with an arbitrary barrier using numerical solutions to Schrodinger's equation.

Main question: what method is best for generating a robust numerical solution to Schrodinger equation for tunneling through an arbitrary barrier? A resource link would be great.

Current solution: I have a time-dependent numerical solution already made using a moving Gaussian wave packet toward a potential barrier.

The tunneling probability is computed as $$T_{prob}=\frac{\Sigma\psi_3^2}{\Sigma\psi_1^2}$$ where $\psi_1, \psi_2, \psi_3$ are the wave function in regions 1, 2, and 3. Region 1 is before the barrier where $V=0$, region 2 is inside the barrier where $E<V$, and region 3 is after the barrier where again $V=0$.

I've discovered, however, that the tunneling probability is dependent on the wavepacket starting position and of course the width $\sigma$ of the Gaussian wave packet. Is there some way to solve this? Wavepacket broadening is expected in a time-dependent solution like this.

I am just looking for a robust method. Perhaps I am missing a concept with this simulation, or there is a better method. The resources I have found almost exclusively focus on bound systems, not tunneling.

Some resources used: https://arxiv.org/html/2403.13857v1#S3 https://www.reed.edu/physics/faculty/wheeler/documents/Quantum%20Mechanics/Miscellaneous%20Essays/Gaussian%20Wavepackets.pdf

Answer 1

I would consider stationary Schroedinger equation instead of solving TDSE: $$ -\frac{1}{2}\frac{\partial^2 \psi}{\partial x^2} = (E - V(x))\psi $$ For large $|x|$, the solution can be written as $$ \psi = \begin{cases} e^{-ikx}, & x < 0\\ ae^{-ikx} + b e^{ikx}, & x > 0,\\ \end{cases} $$ where reflection and transmission amplitudess are $t = 1/a$, $r = b/a$. (It is assumed that the particle is falling on the potential from the right.)

Such a solution can be found by considering the Cauchy problem for the stationary Schroedinger equation with the following initial conditions: $$ \psi(x_0) = 1, $$ $$ \frac{\partial\psi}{\partial x}(x_0) = -ik, $$ where $x_0$ is some large negative coordinate away from the scattering region. It is easy to verify that these conditions indeed correspond to a plane wave $e^{-ikx}$. The coefficients $a$ and $b$ can be extracted from the solution at large positive $x$.

Some drawbacks of this method:

1. The potential should quickly decay away from $x=0$ (ideally, become zero at $|x|$ larger than some $L$)
2. If the potential is high enough, $\psi(x)$ can grow very rapidly in the classically forbidden region, so numerical difficulties may arise.