Table des matières PDFPython

Équation de Schrödinger à une dimension

1. Introduction

Ce document montre comment faire la résolution numérique de l'équation de Schrödinger dépendante du temps, en suivant une méthode similaire à celle utilisée pour l'équation de diffusion. On verra comment appliquer la méthode à un paquet d'onde diffusé par un potentiel.

Voir aussi la simulation Équation de Schrödinger à une dimension.

2. Équation de Schrödinger

On considère l'équation de Schrödinger dépendante du temps pour une fonction d'onde ψ(x,t) :

iψt=-22m2ψx2+V(x)ψ(1)

On s'intéresse à la diffusion de particules dans un potentiel V(x). La résolution est limitée à l'intervalle fini [0,L]. Sur les bords de cet intervalle, la fonction d'onde est supposée nulle. La condition de normalisation s'écrit :

0Lψψ*dx=1(2)

La fonction d'onde initiale ψ(x,0) est connue. Elle prendra la forme d'un paquet d'onde localisée à l'intérieur de l'intervalle [0,L].

On cherche à résoudre numériquement l'équation par la méthode des différences finies. La méthode utilisée est similaire à celle développée dans Équation de diffusion à une dimension.

Soit V0 l'échelle du potentiel. La fonction V(x) utilisée sera sans dimensions. On introduit l'échelle de temps :

τ=V0(3)

et l'échelle de longueur :

λ=2mV0(4)

En utilisant un temps et un x réduits par ces deux échelles, on obtient l'équation :

iψt=-2ψx2+V(x)ψ(5)

On utilise l'opérateur hamiltonien :

H=-2x2+V(x)(6)

La solution de l'équation de Schrödinger peut s'écrire formellement à l'aide de l'exponentielle de l'opérateur hamiltonien ([1]) :

ψ(x,t)=e-iHtψ(x,0)(7)

La conservation de la condition (2) au cours du temps vient du caractère unitaire de l'opérateur e-iHt .

3. Méthode numérique

3.a. Discrétisation

Soit N le nombre de points dans l'intervalle [0,L]. On définit le pas d'espace :

Δx=LN-1(8)

Soit Δt le pas de temps. On pose :

ψjn=ψ(jΔx,nΔt)(9)Vj=V(jΔx)(10)

n est un entier positif ou nul repésentant le temps.

La dérivée temporelle est discrétisée de la manière suivante :

ψtψjn+1-ψjnΔt(11)

La dérivée seconde par rapport à x est discrétisée par :

2ψx2ψj+1n-2ψjn+ψj-1n+1(Δx)2(12)

Le schéma explicite qui en découle est :

iψjn+1-ψjnΔt=-ψj+1n-2ψjn+ψj-1n(Δx)2+Vjψjn(13)

Comme pour l'équation de diffusion, ce schéma explicite a une condition de stabilité très contraignante. On lui préfère donc un schéma implicite.

Pour alléger les notations, on introduit l'opérateur hamiltonien discret défini par :

Hψjn=-ψj+1n-2ψjn+ψj-1n(Δx)2+Vjψjn(14)

3.b. Schéma implicite unitaire

En suivant la démarche adoptée pour l'équation de diffusion à une dimension, on est amené à considérer le schéma implicite suivant :

iψjn+1-ψjnΔt=Hψjn+1(15)

ou encore :

(1+iΔtH)ψjn+1=ψjn(16)

Ce schéma conduit à la résolution, pour chaque pas de temps, d'un système linéaire tridiagonal, pour lequel il existe une méthode directe efficace. Il a toutefois l'inconvénient de ne pas conserver la norme de la fonction d'onde.

Pour un intervalle de temps, on a :

ψ(x,t+Δt)=e-iHΔtψ(x,t)(17)

Le schéma défini ci-dessus revient à utiliser l'approximation au premier ordre suivante :

e-iHΔt1-iHΔt(18)

Cette transformation n'est pas unitaire, c'est pourquoi la norme de la fonction d'onde n'est pas conservée. On la remplace par l'approximation de Cayley ([2]), qui est unitaire :

e-iHΔt1-i12HΔt1+i12HΔt(19)

Le schéma numérique s'écrit alors :

(1+i12ΔtH)ψjn+1=(1-i12ΔtH)ψjn(20)

On obtient ainsi le schéma implicite suivant :

ψj-1n+1+(iα-(Δx)2Vj-2)ψjn+1+ψj+1n+1=-ψj-1n+(iα+(Δx)2Vj+2)ψjn-ψj+1n(21)

avec :

α=2(Δx)2Δt(22)

On le met sous la forme générale suivante :

ajψj-1n+1+bjψjn+1+cjψj+1n+1=djψj-1n+ejψjn+fjψj+1n+sjn(23)

Pour chaque pas de temps, le système tridiagonal est résolu par l'algorithme d'élimination décrit dans Équation de diffusion à une dimension.

Les conditions limites sont de type Dirichlet, puisqu'elles consistent à annuler la fonction d'onde sur les bords de l'intervalle :

ψ0n=ψN-1n=0(24)

3.c. Paquet d'onde

Pour obtenir l'évolution d'un paquet d'onde, on considère comme condition initiale un paquet d'onde gaussien défini par :

ψ(x,0)=eik0xe-(x-x0)24σ02(25)

Le paquet d'onde est initialement en x=x0 et se déplace dans le sens de l'axe x. La densité de propabilité (non normée) associée à ce paquet est :

ρ(x,0)=ψ(x,0)ψ(x,0)*=e-(x-x0)22σ02(26)

La variance de cette densité de probabilité est σ0. Le spectre en nombre d'onde de ce paquet a une largeur donnée par :

σ0Δk=12(27)

Considérons l'évolution du paquet d'onde en l'absence de potentiel (particule libre). La largeur Δk reste constante au cours du temps. Notons σ la largeur du paquet à l'instant t. L'inégalité d'Heisenberg s'écrit :

σΔk12(28)

Pour ce paquet gaussien, l'égalité est réalisée à l'instant t=0. On peut calculer précisément la largeur du paquet à l'instant t ([3]) :

σ=σ01+2t24m2σ04(29)

Cette relation montre l'étalement du paquet d'onde lorsqu'on s'éloigne de l'origine du temps.

La vitesse du paquet d'onde (la vitesse de groupe) est :

v0=k0m(30)

Soit T la durée de la simulation. La distance parcourue par le paquet pendant cette durée doit être de l'ordre de la largeur L. Avec les unités utilisées (m=1/2 et ℏ=1), la condition s'écrit :

TL2k0(31)

Si l'on souhaite que l'étalement du paquet d'onde soit négligeable pendant cette durée, il faut respecter la condition suivante :

Tσ0(32)

qui s'écrit aussi :

2σ02k0L(33)

Par exemple, pour L=1 et σ0=0.05, il faut que k0>200.

Voyons comment choisir l'échantillonnage spatial et temporel. Le pas d'espace doit vérifier :

k0Δx1(34)

La condition de Nyquist-Shannon impose une valeur maximale égale à π, mais il faut prévoir une marge à cause de la dispersion Δk. L'évolution dans le temps de la fonction d'onde se fait à la pulsation moyenne :

ω=E=k022m(35)

Avec les unités utilisées, le pas de temps doit donc vérifier :

k02Δt1(36)

Pour fixer les pas de temps et d'espace, on peut adopter la relation suivante :

Δt=2Δx2(37)

Il suffit alors de vérifier la relation :

k0Δx1(38)

Par exemple, pour k0=200 et L=1, il faut N=2000 points. Le nombre de pas de temps nécessaire est :

P=TΔt=L4k0Δx2(39)

Avec les valeurs précédentes P=5000.

On voit donc que la simulation d'un paquet d'onde présentant un faible étalement est très exigeante en terme de quantité de calculs.

4. Implémentation python

import numpy
import math
import cmath


class Schrodinger1D:
    def __init__(self,L,N):
        self.L = L
        self.N = N
        self.dx = L/(N-1)
        self.b = numpy.zeros(N,dtype=numpy.complex)
        self.e = numpy.zeros(N,dtype=numpy.complex)
        self.x = numpy.zeros(N,dtype=numpy.complex)
        self.beta = numpy.zeros(N,dtype=numpy.complex)
        self.gamma = numpy.zeros(N,dtype=numpy.complex)
        
    def config(self,dt,V):
        N = self.N
        self.dt = dt
        dx2 = self.dx*self.dx
        alpha = 2*dx2/dt
        for k in range(1,self.N-1):
            self.b[k] = 1j*alpha-dx2*V[k]-2.0
            self.e[k] = 1j*alpha+dx2*V[k]+2.0
        self.b[0] = 1.0
        self.b[N-1] = 1.0
        self.beta[0] = self.b[0]
        self.gamma[0] = 1.0/self.beta[0]
        for k in range(1,N):
            self.beta[k] = self.b[k]-self.gamma[k-1]
            if self.beta[k]==0:
                raise Exception("Impossible de resoudre le systeme ligne "+str(k))
            self.gamma[k] = 1.0/self.beta[k]
            
    def iterations(self,psi0,ti,tf):
        if psi0.size!=self.N:
            raise Exception("Taille de U incorrecte")
        psi = psi0.copy()
        t = ti
        while t<tf:
            self.x[0] = (self.e[0]*psi[0]-psi[1])/self.beta[0]
            for k in range(1,self.N-1):
                self.x[k] = (-psi[k-1]+self.e[k]*psi[k]-psi[k+1]-self.x[k-1])/self.beta[k]
            k = self.N-1
            self.x[k] = (-psi[k-1]+self.e[k]*psi[k]-self.x[k-1])/self.beta[k]
            psi[self.N-1] = self.x[self.N-1]
            for k in range(self.N-2,-1,-1):
                psi[k] = self.x[k]-self.gamma[k]*psi[k+1]
            t += self.dt
        return [psi,t]
            
    def V_marche(self,V0):
        V = numpy.zeros(self.N)
        P = int(self.N/2)
        V[P:self.N-1] = V0
        return V
        
    def V_barriere(self,V0,largeur):
        q = int(largeur/self.dx/2)
        V = numpy.zeros(self.N)
        P = int(self.N/2)
        V[P-q:P+q] = V0
        return V
     
    def paquet(self,x0,k0,sigma0):
        psi = numpy.zeros(self.N,dtype=numpy.complex)
        for k in range(1,self.N-1):
            x = self.dx*k
            psi[k] = cmath.exp(1j*k0*x)*math.exp(-(x-x0)*(x-x0)/(4*sigma0**2))
        return psi
            

5. Exemples

5.a. Diffusion de particules sur une marche de potentiel

import numpy
from Schrodinger1D import Schrodinger1D
from matplotlib.pyplot import *

N=2000
L=1.0
dx=L/(N-1)
x = numpy.arange(N)*dx
schrodinger = Schrodinger1D(L,N)
V = schrodinger.V_marche(1.0)
l0 = 2e-2
k0 = numpy.pi*2/l0
E=k0**2
x0 = 0.25
sigma0 = 0.04
psi = schrodinger.paquet(x0,k0,sigma0)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()


                
figAfigA.pdf

On trace le module de la fonction d'onde (densité de probabilité) pour différents instants.

t=0.0
dt=2*dx*dx
v=2*k0
T = 1.0/v
V0 = E
schrodinger.config(dt,V*V0)
[psi,t] = schrodinger.iterations(psi,t,T*0.2)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figBfigB.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.3)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figCfigC.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.5)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figDfigD.pdf

5.b. Diffusion de particules sur une barrière de potentiel

N=2000
L=1.0
dx=L/(N-1)
x = numpy.arange(N)*dx
schrodinger = Schrodinger1D(L,N)
V = schrodinger.V_barriere(1.0,0.05)
l0 = 2e-2
k0 = numpy.pi*2/l0
E=k0*k0
x0 = 0.25
sigma0 = 0.04
psi = schrodinger.paquet(x0,k0,sigma0)
                

On trace le module de la fonction d'onde (densité de probabilité) pour différents instants.

t=0.0
dt=2*dx*dx
v=2*k0
T = 1.0/v
V0=E
schrodinger.config(dt,V*V0)
[psi,t] = schrodinger.iterations(psi,t,T*0.2)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figEfigE.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.3)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figFfigF.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.5)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figGfigG.pdf
[psi,t] = schrodinger.iterations(psi,t,T*0.55)

figure(figsize=(12,6))
plot(x,psi*numpy.conj(psi),'r')
plot(x,V,'b')
xlabel("x")
ylabel("|psi|^2")
grid()
                
figHfigH.pdf
Références
[1]  J.J. Sakurai,  Modern quantum mechanics,  (Addison-Wesley, 1985)
[2]  A. Goldberg, H.M. Schey, J.L. Schwartz,  Computer-generated motion pictures of one-dimensional quantum-mechanical transmission and reflection phenomena,  (Am. J. Phys, 35(3), 1967)
[3]  C. Cohen-Tannoudji, B. Diu, F. Laloë,  Mécanique quantique,  (Hermann, 1980)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.