Table des matières Python

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 :

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

potentiel.svgFigure 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 :

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

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

La solution générale de l'équation est :

avec :

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, kj est imaginaire et il s'agit de solutions exponentielles :

Pour la dernière région BN=0 car on exclue une solution exponentielle croissante et une onde progressant dans le sens de x décroissant. Par ailleurs, la constante A1 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 :

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 :

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

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 :

avec :

Notons M2,1 cette matrice. La matrice inverse M1,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 :

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

qui donne le coefficient de réflexion :

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

La probabilité de transmission est T=1-R. En posant A1=1, le coefficient de réflexion permet de calculer B1. La matrice M21 permet alors de calculer les coefficients (A2,B2), puis la matrice M32 permet de calculer (A3,B3), et ainsi de suite jusqu'aux coefficients (AN,BN). Le coefficient de transmission s'obtient aussi directement par

3. Programme Python

On pose m=1 et . Pour la particule incidente libre (V1=0), on a alors :

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)
        
            

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()
            
figA.svgFigure pleine page

Voici le tracé de la partie réelle de la fonction d'onde pour 2 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()
            
figB.svgFigure 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()
            
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.