Table des matières PDFPython

Particule confinée à une dimension

1. Introduction

Ce document montre comment obtenir les états stationnaires pour une particule confinée dans un potentiel unidimensionnel V(x).

On considère l'équation de Schrödinger pour la fonction d'onde de la particule

On recherche les états stationnaires d'énergie E, définis par :

On doit alors résoudre l'équation de Schrödinger indépendante du temps :

La fonction d'onde recherchée est donc une fonction propre de l'opérateur hamiltonien défini par :

Les énergies E possibles sont des valeurs propres de cet opérateur.

2. Méthode numérique

On pose m=1 et .

On suppose que la particule reste confinée dans un intervalle [a,b]. La fonction d'onde est donc nulle en dehors de cet intervalle, et sur les bornes a et b. On considère un échantillonnage régulier de N points xn sur cet intervalle. On note Δx le pas de x de cette subdivision.

La discrétisation de la dérivée seconde par différences finies donne :

Sur le bord x=a, correspondant à l'indice n=0, on applique la condition limite ϕ-1=0. De même, sur le bord x=b, on applique ϕN=0. On pose :

L'opérateur hamiltonien appliqué à la fonction d'onde discrétisée est :

Il s'agit d'une matrice tridiagonale symétrique.

En notant la matrice colonne contenant les échantillons , on a finalement :

Les énergies sont donc les valeurs propres de la matrice H. Les solutions stationnaires sont les vecteurs propres correspondants.

3. Fonction python

La fonction suivante prend en argument les échantillons de x et du potentiel V, et le nombre de valeurs propres souhaités (les énergies les plus basses). Elle renvoit les énergies et les vecteurs propres correspondants. Les fonctions utilisées sont disponibles avec scipy (version 0.11 ou plus) compilé avec le support de BLAS et LAPACK.

import numpy
import scipy.sparse as sparse
import scipy.sparse.linalg as linalg
from matplotlib.pyplot import *
             
def etats(x,V,nE,tol=0):
    N = V.size
    if x.size!=N:
        raise Exception("tailles de x et V incompatibles")
    alpha = 1.0/(x[1]-x[0])**2
    H = sparse.diags([-alpha/2,V+alpha,-alpha/2],[-1,0,1],shape=(N,N))
    E,phi = linalg.eigs(H,nE,which='SM',tol=tol)
    return numpy.real(E),phi
             

4. Exemples

4.a. Particule dans un puits infini carré

Le potentiel est nul sur tout l'intervalle.

x = numpy.linspace(-1.0,1.0,500)
V = numpy.zeros(x.size)
E,phi = etats(x,V,4)
figure(figsize=(12,8))
for i in range(E.size):
    plot(x,numpy.square(numpy.absolute(phi[:,i])),label="E%d=%f"%(i,E[i]))
xlabel("x")
ylabel("|psi|^2")
grid()
legend(loc='upper right')

                
figA.svgFigure pleine page

On retrouve les états stationnaires pour une particule dans une boîte, identiques aux modes propres d'une corde vibrante (ou d'une cavité).

On trace aussi la répartition des énergies :

nE = 10
E,phi = etats(x,V,nE)
figure()
plot(range(1,nE+1),E,'o')
xlabel("n")
ylabel("E")
grid()
                
figB.svgFigure pleine page

Les énergies sont réparties en n2 :

4.b. Oscillateur harmonique

Le potentiel d'un oscillateur harmonique est :

On se place sur l'intervalle [-1,1]. En x=-1 et x=1 la fonction d'onde est fixée à zéro, ce qui revient à prendre un potentiel infini en ces points. La constante ω doit être choisie assez grande pour que les états obtenus aient leur énergie dans la partie harmonique du potentiel.

omega = 100
x = numpy.linspace(-1.0,1.0,500)
V = x*x*0.5*omega**2
E,phi = etats(x,V,4)
figure(figsize=(12,8))
for i in range(E.size):
    plot(x,numpy.square(numpy.absolute(phi[:,i])),label="E%d=%f"%(i,E[i]))
xlabel("x")
ylabel("|psi|^2")
grid()
legend(loc='upper right')
axis([-0.5,0.5,0,0.025])

                  
figC.svgFigure pleine page
E,phi = etats(x,V,10)
figure()
plot(E,'o')
xlabel("n")
ylabel("E")
grid()
                
figD.svgFigure pleine page

Les niveaux d'énergie sont régulièrement espacés de ω. On retrouve la répartition des énergies de l'oscillateur harmonique :

4.c. Puit de potentiel double

x = numpy.linspace(-1.0,1.0,500)
V = numpy.zeros(x.size)
nL = int(x.size*0.1)
m = int(x.size/2)
V[m-nL:m+nL] = 100.0
figure()
plot(x,V)
xlabel("x")
ylabel("V")
grid()
                
figE.svgFigure pleine page
E,phi = etats(x,V,15)
                

On trace la partie réelle de la fonction d'onde et son module au carré pour les deux premiers états :

figure(figsize=(12,8))
i=0
plot(x,numpy.real(phi[:,i]),label="E%d=%f"%(i,E[i]))
i=1
plot(x,numpy.real(phi[:,i]),label="E%d=%f"%(i,E[i]))
xlabel("x")
ylabel("Re(psi)")
grid()
legend(loc='upper center')
                
figF.svgFigure pleine page
figure(figsize=(12,8))
i=0
plot(x,numpy.square(numpy.absolute(phi[:,i])),label="E%d=%f"%(i,E[i]))
i=1
plot(x,numpy.square(numpy.absolute(phi[:,i])),label="E%d=%f"%(i,E[i]))
xlabel("x")
ylabel("|psi|^2")
grid()
legend(loc='upper center')
                
figG.svgFigure pleine page

Les deux premiers états des énergies très proches. Pour le premier, la fonction d'onde est symétrique. Pour le deuxième, elle est antisymétrique. La densité de probabilité est la même pour les deux états, avec une particule délocalisée sur les deux puits.

En faisant la somme (ou la différence) de ces deux états, on obtient une particule localisée seulement dans un des puits (état non stationnaire) :

phiD = phi[:,0]+phi[:,1]
phiG = phi[:,0]-phi[:,1]
figure()
plot(x,numpy.square(numpy.absolute(phiD)),label="psi0+psi1")
plot(x,numpy.square(numpy.absolute(phiG)),label="psi0-psi1")
xlabel("x")
ylabel("|psi|^2")
grid()
legend(loc='upper center')
                  
figH.svgFigure pleine page

Voici les deux états suivants :

figure(figsize=(12,8))
i=2
plot(x,numpy.real(phi[:,i]),label="E%d=%f"%(i,E[i]))
i=3
plot(x,numpy.real(phi[:,i]),label="E%d=%f"%(i,E[i]))
xlabel("x")
ylabel("Re(psi)")
grid()
legend(loc='upper center')
                
figI.svgFigure pleine page

Ces deux niveaux n'apparaissent pas dans l'ordre d'énergie croissante. L'état symétrique est toujours celui de plus basse énergie.

figure(figsize=(12,8))
i=2
plot(x,numpy.square(numpy.absolute(phi[:,i])),label="E%d=%f"%(i,E[i]))
i=3
plot(x,numpy.square(numpy.absolute(phi[:,i])),label="E%d=%f"%(i,E[i]))
xlabel("x")
ylabel("|psi|^2")
grid()
legend(loc='upper center')
                
figJ.svgFigure pleine page

Voici les niveaux d'énergie :

E = sorted(E)
figure()
plot(E,'o')
xlabel("n")
ylabel("E")
grid()
                
figK.svgFigure pleine page

Seuls les 8 premiers niveaux ont une énergie inférieure à la hauteur de la barrière séparant les deux puits. Pour les énergies supérieures à la hauteur, on retrouve la répartition d'un puit simple. En augmentant la hauteur de la barrière, on peut bien sûr augmenter le nombre d'états d'énergie inférieure à cette hauteur.

4.d. Potentiel périodique

On modifie la matrice H pour appliquer une condition limite périodique sur la fonction d'onde.

def etats_periodique(x,V,nE,tol=0):
    N = V.size
    if x.size!=N:
        raise Exception("tailles de x et V incompatibles")
    alpha = 1.0/(x[1]-x[0])**2
    H = sparse.diags([-alpha/2,V+alpha,-alpha/2],[-1,0,1],shape=(N,N))
    H = H.toarray()
    H[0,N-1] = -alpha/2
    H[N-1,0] = -alpha/2
    E,phi = linalg.eigs(H,nE,which='SM',tol=tol)
    return numpy.real(E),phi
             

Voici un exemple avec un potentiel sinusoïdal :

x = numpy.linspace(-1.0,1.0,500)
V = (numpy.sin(40*numpy.pi*x)+1)*0.5*10000
figure()
plot(x,V)
xlabel("x")
ylabel("V")
grid()
                
figL.svgFigure pleine page
E,phi = etats_periodique(x,V,100)
E = sorted(E)
figure()
plot(E,'o')
xlabel("n")
ylabel("E")
grid()
                
figM.svgFigure pleine page

Les énergies se répartissent par bandes. Les bandes d'énergie s'observent pour les électrons de conduction dans un métal. Ces électrons sont en effet soumis à un potentiel périodique de la part des cations du réseau cristallin.

figure(figsize=(12,8))
i=0
plot(x,numpy.real(phi[:,i]),label="E%d=%f"%(i,E[i]))
i=1
plot(x,numpy.real(phi[:,i]),label="E%d=%f"%(i,E[i]))
xlabel("x")
ylabel("Re(psi)")
grid()
legend(loc='upper center')
                
figN.svgFigure pleine page
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.