Diffusion d'une particule dans un potentiel constant par morceaux

1. Problème

Une particule matérielle (non relativiste) de masse m se déplace dans un espace unidimensionnel. On considère l'équation de Schrödinger pour sa fonction d'onde $ψ (x, t)$ :

i ℏ \frac{\partial ψ}{\partial t} = - \frac{ℏ^{2}}{2 m} \frac{\partial^{2} ψ}{\partial x^{2}} + V (x) ψ (1)

On suppose que la particule se déplace initialement dans le sens de x croissant sans potentiel (V=0). Elle rencontre, à partir de x₁ une région où le potentiel est constant par morceaux.

Figure pleine page

Il y a N régions et N-1 discontinuités du potentiel.

On s'intéresse à des états stationnaires de la particule, d'énergie E :

ψ (x, t) = e^{- i \frac{E}{ℏ} t} φ (x) (2)

On est donc amené à résoudre l'équation de Schrödinger indépendante du temps :

\frac{d^{2} φ_{j} (x)}{d x^{2}} + \frac{2 m (E - V_{j})}{ℏ^{2}} φ_{j} (x) = 0 (3)

pour chaque intervalle [x_j-1,x_j], avec les conditions de continuité aux frontières :

\begin{matrix} φ_{k} (x_{j}) = φ_{j + 1} (x_{j}) & (4) \\ \frac{d φ_{j}}{d x} (x_{j}) = \frac{d φ_{j + 1}}{d x} (x_{j}) & (5) \end{matrix}

La solution générale de l'équation $(3)$ est :

φ_{j} (x) = A_{j} e^{i k_{j} x} + B_{j} e^{- i k_{j} x} (6)

avec :

ℏ k_{j} = \sqrt{2 m (E - V_{j})} (7)

Si l'énergie est supérieure au potentiel, il s'agit d'une solution sinusoïdale, somme de deux ondes progressives de sens opposés. Si l'énergie est inférieure au potentiel, k_j est imaginaire et il s'agit de solutions exponentielles :

\begin{matrix} k_{j} = i α_{j} & (8) \\ φ_{j} (x) = A_{j} e^{- α_{j} x} + B_{j} e^{α_{j} x} & (9) \end{matrix}

Pour la dernière région B_N=0 car on exclue une solution exponentielle croissante et une onde progressant dans le sens de x décroissant. Par ailleurs, la constante A₁ pour l'onde incidente est arbitraire. Il reste donc 2(N-2)+2=2N-2 constantes à déterminer. Les conditions aux frontières fournissent 2(N-1)=2N-2 équations linéaires pour ces constantes, ce qui permet de les déterminer toutes.

2. Formulation matricielle

Considérons par exemple la frontière entre les régions 1 et 2. Les fonctions d'onde sont :

\begin{matrix} φ_{1} (x) = A_{1} e^{i k_{1} x} + B_{1} e^{- i k_{1} x} & (10) \\ φ_{2} (x) = A_{2} e^{i k_{2} x} + B_{2} e^{- i k_{2} x} & (11) \end{matrix}

On choisit de faire les calculs avec le nombre d'onde k, qui est imaginaire pur lorsque l'énergie est inférieure au potentiel.

Pour alléger les notations, on introduit :

\begin{matrix} e_{1} = e^{i k_{1} x_{1}} & (12) \\ f_{1} = e^{- i k_{1} x_{1}} & (13) \\ e_{2} = e^{i k_{2} x_{1}} & (14) \\ f_{2} = e^{- i k_{2} x_{1}} & (15) \end{matrix}

La continuité de la fonction d'onde et de sa dérivée s'écrivent :

\begin{matrix} A_{1} e_{1} + B_{1} f_{1} = A_{2} e_{2} + B_{2} f_{2} & (16) \\ k_{1} (A_{1} e_{1} - B_{1} f_{1}) = k_{2} (A_{2} e_{2} - B_{2} f_{2}) & (17) \end{matrix}

La résolution de ce système donne les constantes de la région 2 en fonction de celle de la région 1, sous forme matricielle :

(\begin{matrix} A_{2} \\ B_{2} \end{matrix}) = (\begin{matrix} a & b \\ c & d \end{matrix}) (\begin{matrix} A_{1} \\ B_{1} \end{matrix}) (18)

avec :

\begin{matrix} a = \frac{(k_{1} + k_{2}) e_{1}}{2 k_{2} e_{2}} & (19) \\ b = \frac{(k_{2} - k_{1}) f_{1}}{2 k_{2} e_{2}} & (20) \\ c = \frac{(k_{2} - k_{1}) e_{1}}{2 k_{2} f_{2}} & (21) \\ d = \frac{(k_{1} + k_{2}) f_{1}}{2 k_{2} f_{2}} & (22) \end{matrix}

Notons M_2,1 cette matrice. La matrice inverse M_1,2, qui permet d'exprimer les constantes de la région 1 en fonction de celles de la région 2 se déduit de la première par une permutation des indices 1 et 2.

Pour le passage de la région 1 à la dernière, on écrit un produit matriciel :

(\begin{matrix} A_{N} \\ B_{N} \end{matrix}) = M_{N, N - 1} M_{N - 1, N - 2} . . . M_{21} (\begin{matrix} A_{1} \\ B_{1} \end{matrix}) = (\begin{matrix} a & b \\ c & d \end{matrix}) (\begin{matrix} A_{1} \\ B_{1} \end{matrix}) (23)

La constante B_N étant nulle, on a l'équation :

0 = c A_{1} + d B_{1} (24)

qui donne le coefficient de réflexion :

r = \frac{B_{1}}{A_{1}} = - \frac{c}{d} (25)

d'où l'on déduit la probabilité de réflexion :

R = | r |^{2} (26)

La probabilité de transmission est T=1-R. En posant A₁=1, le coefficient de réflexion permet de calculer B₁. La matrice M₂₁ permet alors de calculer les coefficients (A₂,B₂), puis la matrice M₃₂ permet de calculer (A₃,B₃), et ainsi de suite jusqu'aux coefficients (A_N,B_N). Le coefficient de transmission s'obtient aussi directement par

τ = \frac{A_{N}}{A_{1}} = a + b r (27)

3. Programme Python

Pour le traitement numérique, il est préférable d'utiliser des grandeurs sans dimensions. Introduisons la longueur définie par :

ℓ = \frac{ℏ}{\sqrt{2 m E_{0}}} (28)

où E₀ est une énergie de référence. Si l'on pose $x' = x / ℓ$ , E'=E/E₀ et V'=V/E₀, l'équation de Schrödinger s'écrit :

\frac{d^{2} φ_{j} (x')}{d x'^{2}} + (E' - V'_{j}) φ_{j} (x') = 0 (29)

Formellement, cela revient à poser m=1 et $ℏ = 1$ dans l'équation initiale.

Si la particule est un électron si l'on choisit E₀=1 eV on a $ℓ = 0, 19 n m$ .

DiffusionParticule.py

import numpy
import cmath
 
class DiffusionParticule:
    def __init__(self,x,V):
        self.N = len(V)
        self.V = V
        self.x = x
        if len(x)!=self.N-1:
            raise Exception("Tailles de x et V icompatibles")
        
        
    
    def reflexion(self,E):
        self.E = E
        self.M = []
        self.k = []
        for j in range(self.N-1):
            k1 = cmath.sqrt(2*(E-self.V[j]))
            k2 = cmath.sqrt(2*(E-self.V[j+1]))
            self.k.append(k1)
            e1 = cmath.exp(1j*k1*self.x[j])
            f1 = cmath.exp(-1j*k1*self.x[j])
            e2 = cmath.exp(1j*k2*self.x[j])
            f2 = cmath.exp(-1j*k2*self.x[j])
            a = (k1+k2)*e1/(2*k2*e2)
            b = (k2-k1)*f1/(2*k2*e2)
            c = (k2-k1)*e1/(2*k2*f2)
            d = (k1+k2)*f1/(2*k2*f2)
            self.M.append(numpy.array([[a,b],[c,d]],dtype=numpy.complex))
        self.k.append(k2)
        self.M_tot = self.M[0]
        for j in range(1,self.N-1):
            self.M_tot = self.M[j].dot(self.M_tot)
        a = self.M_tot[0,0]
        b = self.M_tot[0,1]
        c = self.M_tot[1,0]
        d = self.M_tot[1,1]
        r = -c/d
        t = a+b*r
        R = abs(r)**2
        A = 1
        B = r
        self.A = [A]
        self.B = [B]
        for j in range(self.N-1):
            AB = self.M[j].dot(numpy.array([A,B]))
            A = AB[0]
            B = AB[1]
            self.A.append(A)
            self.B.append(B)
        self.B[self.N-1] = 0
        return R
        
    def echantillons(self,xmin,xmax,t,ne):
        dx = (xmax-xmin)*1.0/ne
        x = numpy.zeros(0)
        psi = numpy.zeros(0)
        x1 = xmin
        for j in range(self.N-1):
            x2 = self.x[j]
            n = (x2-x1)/dx
            x12 = numpy.linspace(x1,x2,n)
            x1 = x2
            psi = numpy.append(psi,self.A[j]*numpy.exp(1j*self.k[j]*x12)+self.B[j]*numpy.exp(-1j*self.k[j]*x12))
            x = numpy.append(x,x12)
        x2 = xmax
        n = (x2-x1)/dx
        x12 = numpy.linspace(x1,x2,n)
        j = self.N-1
        psi = numpy.append(psi,self.A[j]*numpy.exp(1j*self.k[j]*x12)+self.B[j]*numpy.exp(-1j*self.k[j]*x12))
        x = numpy.append(x,x12)
        psi = psi*numpy.exp(-1j*self.E*t)
        return (x,psi)

4. Exemples

4.a. Barrière de potentiel carrée

Voici par exemple une barrière de potentiel et le tracé du coefficient de transmission en fonction de l'épaisseur de la barrière.

from DiffusionParticule import *
from matplotlib.pyplot import *

L = numpy.linspace(0,5.0,100)
T0 = numpy.zeros(L.size)
T1 = numpy.zeros(L.size)
T2 = numpy.zeros(L.size)
E0 = 0.1
E1 = 0.5
E2 = 0.8
for i in range(L.size):
    diffusion = DiffusionParticule([0,L[i]],[0,1,0])
    T0[i] = 1-diffusion.reflexion(E0)
    T1[i] = 1-diffusion.reflexion(E1)
    T2[i] = 1-diffusion.reflexion(E2)
    
figure()
plot(L,T0,label="E=0.1")
plot(L,T1,label="E=0.5")
plot(L,T2,label="E=0.8")
xlabel("L")
ylabel("T")
legend(loc = "upper right")
grid()

Figure pleine page

Voici le tracé de la partie réelle de la fonction d'onde pour deux instants :

diffusion = DiffusionParticule([0,1],[0,1,0])
E = 0.5
diffusion.reflexion(E)
xmin=-20
xmax=20

figure()
(x,psi) = diffusion.echantillons(xmin,xmax,0.0,1000)
plot(x,numpy.real(psi))
(x,psi) = diffusion.echantillons(xmin,xmax,1.0,1000)
plot(x,numpy.real(psi))
(x,psi) = diffusion.echantillons(xmin,xmax,2.0,1000)
plot(x,numpy.real(psi))
xlabel("x")
ylabel("Re(psi)")
grid()

Figure pleine page

Le code suivant permet de faire une animation :

import matplotlib.animation as animation

(x,psi) = diffusion.echantillons(xmin,xmax,0.0,1000)

fig, ax = subplots()
line, = ax.plot(x,numpy.real(psi))
ax.grid()
ax.axis([xmin,xmax,-2,2])
temps = 0.0
dt = 4.0e-2

def animate(i):
    global temps
    temps += dt
    (x,psi) = diffusion.echantillons(xmin,xmax,temps,1000)
    line.set_xdata(x)
    line.set_ydata(numpy.real(psi))
    return line,

ani = animation.FuncAnimation(fig,animate,1000,interval=40)
    
show()

4.b. Barrière de potentiel en cosinus

On considère une barrière de potentiel formée par un arc de sinusoïde (une période). La courbe est échantillonnée pour être représentée par une fonction constante par morceaux.

x = numpy.linspace(-1,1,100)
V = (numpy.cos(x*numpy.pi)+1.0)*0.5
figure()
plot(x,V)
xlabel("x")
ylabel("V")
grid()

Figure pleine page

La largeur à mi-hauteur est L=1, ce qui permet d'obtenir un effet tunnel très important.

On trace le coefficient de transmission en fonction de l'énergie :

x = numpy.delete(x,x.size-1)
diffusion = DiffusionParticule(x,V)
                
E = numpy.linspace(0.01,3.0,100)
T = numpy.zeros(E.size)
for i in range(E.size):
    R = diffusion.reflexion(E[i])
    T[i] = 1-R
    
figure()
plot(E,T)
xlabel("E/V0")
ylabel("T")
grid()

Figure pleine page