Meaning of transmission coefficient greater than one in a potential well problem

Consider the finite 1D wedge-shaped potential well given by

$$V(x)=V_0\left(\frac{|x|}{a}-1\right) \hspace{10pt}\mathrm{for}\hspace{3pt} |x|<a;\hspace{6pt}V(x)=0 \hspace{10pt}\mathrm{for}\hspace{3pt} |x|>a.$$ I'm trying to find the reflection and transmission coefficients for a stream of electrons coming from $x=-\infty$ (so $E>0$, no bound states involved). In order to do that, I've divided the domain into 4 parts (like in the picture), solved the Schrödinger equation on each one, got a linear system for the coefficients of each wavefunction and solved it. Meanwhile I've also solved the problem numerically, and got pretty expected results, like the one below. All was well untill I plugged the analytic solution for the transmission coefficient into Python to render a graph and got this. $E$ is in eV, $x$ and $a$ in nm. Now, I've gone over my calculations a dozen times, and I still can't spot a mistake, so I want to know wheather it's possible to have a transmission coefficient greater than one for certain values of energy, like in this scenario.

I'm aware that this breaks conservation of energy (red alert!), and I'm still hoping that this is due to a numerical error, but I want to hear some other ideas regarding this phenomenon. I've found this as a related question, but I'm not really satisfied with the answers. Why doesn't this phenomenon happen in a rectangular well? Should I discard the $T>1$ solutions and leave energy gaps in the graph? That seems pretty arbitrary arbitrary to me, but it may be the case here.

EDIT: The comments suggest to give an insight into my analytical solution.

Some initial substitutions: $k=\sqrt{\frac{2mE}{\hbar^2}}$, $k_0=\sqrt{\frac{2mV_0}{\hbar^2}}$, $q_0=(k_0a)^{2/3}$, $\epsilon^2=\frac{E}{V_0}$.

The time-independant Schrödinger equation for the wavefunction on each section reads

$$\frac{d^2\psi_I}{dx^2}+k^2\psi_I=0$$ $$\frac{d^2\psi_{II}}{dx^2}+(-\frac{x}{a}+\epsilon^2-1)k_0^2\psi_{II}=0$$ $$\frac{d^2\psi_{III}}{dx^2}+(\frac{x}{a}+\epsilon^2-1)k_0^2\psi_{III}=0$$ $$\frac{d^2\psi_{IV}}{dx^2}+k^2\psi_{IV}=0$$

Solutions to $I$ and $IV$ may be immediatly given as

$$\psi_I(x)=Ae^{ikx}+Be^{-ikx}$$ $$\psi_{IV}(x)=Ge^{ikx},$$

the first corresponding to the incoming and reflected wave, and the other to the transmitted wave.

Introducing the change of variables

$$\xi(x)=q_0(\frac{|x|}{a}+\epsilon^2-1)$$ transforms the other two equations into Airy differential equations:

$$\frac{d^2\psi_{II}}{d\xi^2}=-\xi\psi_{II}$$ $$\frac{d^2\psi_{III}}{d\xi^2}=-\xi\psi_{III}$$

with general solutions in the form

$$\psi_{II}(x)=C\mathrm{Ai}(-\xi)+D\mathrm{Bi}(-\xi)$$ $$\psi_{III}(x)=E\mathrm{Ai}(-\xi)+F\mathrm{Bi}(-\xi),$$

where $\mathrm{Ai}(z)$ and $\mathrm{Bi}(z)$ are the standard Airy functions of the first and second kind.

Using the fact that both $\psi(x)$ and $\psi'(x)$ are continuous, this gives the linear system:

$$\left(\begin{matrix} e^{-ika} & e^{ika} & -\mathrm{Ai}(-\xi_a) & -\mathrm{Bi}(-\xi_a) & 0 & 0 \\ -\frac{ika}{q_0}e^{-ika} & \frac{ika}{q_0}e^{ika} & \mathrm{Ai}'(-\xi_a) & \mathrm{Bi}'(-\xi_a) & 0 & 0 \\ 0 & 0 & \mathrm{Ai}(-\xi_0) & \mathrm{Bi}(-\xi_0) & -\mathrm{Ai}(-\xi_0) & -\mathrm{Bi}(-\xi_0) \\ 0 & 0 & \mathrm{Ai}'(-\xi_0) & \mathrm{Bi}'(-\xi_0) & \mathrm{Ai}'(-\xi_0) & \mathrm{Bi}'(-\xi_0) \\ 0 & 0 & 0 & 0 & \mathrm{Ai}(-\xi_a) & \mathrm{Bi}(-\xi_a) \\ 0 & 0 & 0 & 0 & \mathrm{Ai}'(-\xi_a) & \mathrm{Bi}(-\xi_a) \end{matrix} \right)\left(\begin{matrix} A \\ B \\ C \\ D \\ E \\ F \end{matrix} \right)=\left( \begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ Ge^{ika} \\ -G\frac{ika}{q_0}e^{ika} \end{matrix} \right)$$

where $\xi_a\equiv\xi(x=a)=\xi(-a)=q_0\epsilon^2$ and $\xi_0\equiv\xi(0)=q_0(\epsilon^2-1)$ appear as shorthand notations. Prime corresponds to $\frac{d}{dx}$. The currents are given by

$$J_I=\frac{\hbar}{m}\Im\left(\psi_I^*\frac{d\psi_I}{dx}\right)=\frac{\hbar}{m}(|A|^2-|B|^2)=J_{in}-J_{ref}$$ $$J_{IV}=\frac{\hbar}{m}\Im\left(\psi_{IV}^*\frac{d\psi_{IV}}{dx}\right)=\frac{\hbar}{m}|G|^2=J_{tr}$$

Since $T+R=1$, it follows that $|A|^2-|B|^2=|G|^2$. $G$ is actually undetermined and is used as a free coefficient (every other coefficient can be expressed as $L=G\cdot blabla$), so it may be freely set to be 1. It follows that $T=\frac{J_{tr}}{J_{in}}=\frac{1}{|A|^2}$ and $R=1-\frac{1}{|A|^2}$. Hence, only $A$ is needed to calculate the transmission and reflexion coefficients. It may be calculated from the above system using Cramer's rule.

My result:

\begin{equation*}\label{parametara} \begin{split} A&=\frac{iq_0\pi^2e^{2ika}}{ka} \left[\left(\left(\frac{ka}{q_0}\right)^2\mathrm{Ai}^2(-\xi_a)+\mathrm{Ai}'^2(-\xi_a)\right)\mathrm{Bi}(-\xi_0)\mathrm{Bi}'(-\xi_0)\right.\\ &+\left(\left(\frac{ka}{q_0}\right)^2\mathrm{Bi}^2(-\xi_a)+\mathrm{Bi}'^2(-\xi_a)\right)\mathrm{Ai}(-\xi_0)\mathrm{Ai}'(-\xi_0)\\ &\left.-\left(\left(\frac{ka}{q_0}\right)^2\mathrm{Ai}(-\xi_a)\mathrm{Bi}(-\xi_a)+\mathrm{Ai}'(-\xi_a)\mathrm{Bi}'(-\xi_a)\right)(\mathrm{Ai}(-\xi_0)\mathrm{Bi}'(-\xi_0)+\mathrm{Ai}'(-\xi_0)\mathrm{Bi}(-\xi_0)) \right], \end{split} \end{equation*}

where use was made of the Wronskian for Airy functions:

$$\mathcal{W}[\mathrm{Ai}(z),\mathrm{Bi}(z)]=\pi^{-1}.$$

The next step was to plug the absolute square of $A$ into expressions for $T$ and $R$, which rendered the problematic graph.

• I would suspect that your analytical solution does not satisfy the boundary conditions. So maybe you should give the details of your analytical solution, otherwise it is difficult to say anything for sure. – akhmeteli Dec 3 '16 at 2:28
• You have to specify how you arrive at this T>1 result and give some pertinent formulas. Otherwise it is difficult to say anything meaningful to your problem. – freecharly Dec 3 '16 at 2:53
• Your $T$ result seems pretty close to the expected. Numerical error? – Gert Dec 3 '16 at 3:43
• Edited the question accordingly. – Soba noodles Dec 3 '16 at 12:52
• Would Computational Science be a better home for this question? – Qmechanic Dec 3 '16 at 14:49

I haven't studied your solution in detail, but, for what it's worth, the formulas for the solutions of Airy equation differ depending on the sign of the coefficient. Have you checked that your solutions are applicable independent of the value of $E$? You may wish to check if your solutions give correct results for $V_0=0$.

• I've tried $V_0=0$, but the algorhitm crashes because somewhere a division with zero occurs (since $A\propto V_0$). I've played with other values of $a$ and $V_0$ as well, and the issue of $T>1$ persists, if that's what you mean. – Soba noodles Dec 3 '16 at 15:01
• @Sobanoodles: That's not what I mean. The solution of the Airy equation bifurcates depending on whether you have $-k^2$ or $+k^2$ (in notation of mathworld.wolfram.com/AiryDifferentialEquation.html). As for the crashing algorithm, you can try small $V_0$. – akhmeteli Dec 3 '16 at 15:14
• @akhmetali I see what you mean. Since I've got the Airy equation in the form $\frac{d^2\psi}{d\xi^2}+\xi\psi=0$, I've simply declared that the solutions are of the form $\psi(\xi)=K\mathrm{Ai}(-\xi)+J\mathrm{Bi}(-\xi),$ which is justified as a change of variables $\xi=-z$ in the conventional Airy equation $\frac{d^2\psi}{dz^2}-z\psi=0$ – Soba noodles Dec 3 '16 at 15:21
• As for the small $V_0$ case, here it is. As expected, when there is no potential well, the whole wave gets transmitted. – Soba noodles Dec 3 '16 at 15:59
• @Sobanoodles: Glad I could help. – akhmeteli Dec 4 '16 at 18:01

Change $$\psi_{IV}$$ to $$G e^{-i k x}$$ Calculate the reflection probability instead then find T = 1-R that should work well.

Didn't you forget to set a condition which reflects the fact that there is no incoming particles from the right, that is, the amplitude of the phase $e^{-ixk}$ is zero for $x\to +\infty$?

• I've edited the question now to show how I got to the graph, so you can see that I've set $\psi_{IV}=e^{ikx}$. In other words, I didn't forget about setting that amplitude to 0, but thank you for pointing that out. – Soba noodles Dec 3 '16 at 13:02
• Somehow you managed to transform this problem from triangle well to triangle barrier: – Пегасище Dec 4 '16 at 0:17
• Would you please elaborate that? Has it something to do with the change of variables discussed in akhmeteli's answer and comments? – Soba noodles Dec 4 '16 at 1:15
• Sorry, i pressed enter and the page sent my brocken comment. About trianlge barrier: in the well: $k^2=\frac{2m(E-V)}{\hbar^2} = \frac{2mV_0}{\hbar^2}(\frac{E}{V_0}-\frac{|x|}{a}+1)=k_0^2(\epsilon^2-\frac{|x|}{a}+1)$ so you must rewrite equations for $\psi_{II}$ and $\psi_{III}$. I am completely not sure if this is significant for your problem. – Пегасище Dec 4 '16 at 1:21
• What confuses me more is the sings of the elements in your matrix. Why did you set (1,3), (1,4), (3,5) and (3,6) with -1 multiplier? Also did you mention that for one of $\psi_{II}$ and $\psi_{III}$ : $\frac{d}{dx} = \frac{d}{d\xi}$ and for the other $\frac{d}{dx} = -\frac{d}{d\xi}$ – Пегасище Dec 4 '16 at 1:24