Meccanica quantistica - Dipartimento di Fisica

Appunti per il corso di

Meccanica quantistica

Corso di Laurea in Fisica Computazionale Universit` a di Udine Anno accademico 2009/2010

Paolo Giannozzi1 sulla base del software e delle note scritte da

Furio Ercolessi1 e Stefano de Gironcoli2 1 Universit` a

di Udine - Dipartimento di Fisica 2 SISSA - Trieste

Versione del: October 28, 2010

Contents Prefazione

1

1 Meccanica classica 1.1 Formulazione Lagrangiana e Hamiltoniana della meccanica 1.2 Un’applicazione classica: moto di un punto materiale . . . . 1.2.1 L’algoritmo di Størmer-Verlet . . . . . . . . . . . . . 1.3 Programma: newton . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Struttura del programma . . . . . . . . . . . . . . . 1.3.2 Laboratorio . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

3 3 5 6 6 7 8

2 Introduzione all’equazione di Schr¨ odinger 2.1 Verso la Meccanica Quantistica . . . . . . . . 2.1.1 Dualismo onda-particella: fotoni . . . 2.1.2 Quantizzazione dei livelli di energia . . 2.1.3 Dualismo onda-particella: elettroni . . 2.1.4 Principio di indeterminazione . . . . . 2.2 L’equazione di Schr¨ odinger per una particella 2.2.1 L’equazione di Schrödinger dipendente 2.3 La particella libera . . . . . . . . . . . . . . . 2.4 Pacchetti d’onda . . . . . . . . . . . . . . . . 2.5 Potenziali modello . . . . . . . . . . . . . . . 2.5.1 Gradino di potentiale . . . . . . . . . 2.5.2 Barriera di potenziale . . . . . . . . . 2.5.3 Buca di potenziale . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

10 10 10 11 11 12 12 13 15 16 18 18 19 20

. . . . . . . . . . . . dal . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . tempo . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

3 L’equazione di Schr¨ odinger unidimensionale: soluzione analitica e numerica 3.1 L’oscillatore armonico . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Unit` a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Soluzione e livelli energetici . . . . . . . . . . . . . . . . . 3.1.3 Energia di punto zero . . . . . . . . . . . . . . . . . . . . 3.1.4 Simmetria e parità . . . . . . . . . . . . . . . . . . . . . . 3.1.5 Confronto con la densità di probabilità classica . . . . . . 3.2 Meccanica quantistica e codici numerici: alcune considerazioni . 3.2.1 Energie cinetiche negative . . . . . . . . . . . . . . . . . . 3.2.2 Effetti della quantizzazione . . . . . . . . . . . . . . . . .

i

23 23 24 24 26 26 27 28 28 28

3.3

Il metodo di Numerov . . . . . 3.3.1 Programma: harmonic0 3.3.2 Programma: harmonic1 3.3.3 Laboratorio . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

29 31 32 33

4 Propriet` a dell’equazione di Schr¨ odinger 4.1 Ortonormalit` a delle funzioni d’onda . . 4.2 Sviluppo di una soluzione generica . . . 4.3 Valori medi . . . . . . . . . . . . . . . . 4.4 La formulazione matriciale . . . . . . . . 4.5 Regole di commutazione . . . . . . . . . 4.6 Quantit` a conservate . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

35 35 37 38 39 41 42

5 Atomi con un elettrone 5.1 Equazione di Schr¨ odinger in un campo centrale . 5.2 Il momento angolare . . . . . . . . . . . . . . . . 5.3 Autofunzioni del momento angolare . . . . . . . . 5.4 Separazione in parte radiale e angolare . . . . . . 5.4.1 Funzioni d’onda angolari . . . . . . . . . . 5.5 Il potenziale coulombiano . . . . . . . . . . . . . 5.6 La funzione d’onda radiale per atomi idrogenoidi 5.6.1 Densit` a radiale . . . . . . . . . . . . . . . 5.6.2 Stato fondamentale . . . . . . . . . . . . . 5.6.3 Comportamento vicino al nucleo . . . . . 5.6.4 Comportamento lontano dal nucleo . . . . 5.6.5 Numero di nodi . . . . . . . . . . . . . . . 5.7 Degenerazione accidentale e simmetria dinamica 5.8 Programma: hydrogen . . . . . . . . . . . . . . . 5.8.1 Griglia logaritmica . . . . . . . . . . . . . 5.8.2 Applicazione della teoria perturbativa . . 5.8.3 Laboratorio . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

45 45 46 47 49 51 51 52 54 54 54 55 55 55 56 56 57 59

6 Metodi approssimati 6.1 Metodo perturbativo . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1 Perturbazioni con autovalori degeneri . . . . . . . . . . . 6.2 Perturbazioni dipendenti dal tempo: transizioni elettromagnetiche 6.2.1 Transizioni di dipolo . . . . . . . . . . . . . . . . . . . . . 6.3 Metodo variazionale . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.1 Dimostrazione del principio variazionale (I) . . . . . . . . 6.3.2 Dimostrazione del principio variazionale (II) . . . . . . . . 6.3.3 Energia dello stato fondamentale . . . . . . . . . . . . . . 6.3.4 Il metodo variazionale in pratica . . . . . . . . . . . . . . 6.4 Problema secolare . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4.1 Sviluppo in funzioni ortonormali . . . . . . . . . . . . . . 6.4.2 Sviluppo in funzioni non ortonormali . . . . . . . . . . . . 6.5 Programma: hydrogen gauss . . . . . . . . . . . . . . . . . . . . 6.5.1 Laboratorio . . . . . . . . . . . . . . . . . . . . . . . . . .

ii

60 60 62 62 64 65 66 67 67 68 68 69 71 73 74

6.6 6.7

Base di onde piane . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Programma: pwell . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.7.1 Laboratorio . . . . . . . . . . . . . . . . . . . . . . . . . . 76

7 Atomi a pi` u elettroni 7.1 Lo spin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Composizione di momenti angolari: la rappresentazione accoppiata 7.2.1 Esempio: singoletti e tripletti . . . . . . . . . . . . . . . . 7.2.2 Presenza di accoppiamento . . . . . . . . . . . . . . . . . 7.3 Particelle identiche: principio di indistinguibilità . . . . . . . . . 7.4 Operatori di permutazione . . . . . . . . . . . . . . . . . . . . . . 7.5 Caso di pi` u particelle e sistemi composti . . . . . . . . . . . . . . 7.6 Determinanti di Slater . . . . . . . . . . . . . . . . . . . . . . . . 7.7 Atomi a due elettroni . . . . . . . . . . . . . . . . . . . . . . . . 7.8 Trattamento perturbativo dell’atomo di elio . . . . . . . . . . . . 7.9 Trattamento variazionale dell’atomo di elio . . . . . . . . . . . . 7.10 Programma: helium gauss . . . . . . . . . . . . . . . . . . . . . . 7.10.1 Laboratorio . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Metodo di Hartree-Fock 8.1 Il metodo di Hartree . . . . . . . . . . . . . 8.1.1 Definizioni . . . . . . . . . . . . . . . 8.1.2 Equazioni di Hartree . . . . . . . . . 8.1.3 Significato del potenziale di Hartree 8.1.4 Campo autoconsistente . . . . . . . 8.1.5 Autovalori ed Energia di Hartree . . 8.2 Il metodo di Hartree-Fock . . . . . . . . . . 8.2.1 Potenziale colombiano e di scambio . 8.2.2 La densit` a di scambio . . . . . . . . 8.2.3 L’atomo di elio . . . . . . . . . . . . 8.3 L’energia di correlazione . . . . . . . . . . . 8.4 Programma: helium hf radial . . . . . . . . 8.4.1 Laboratorio . . . . . . . . . . . . . . 8.5 Programma: helium hf gauss . . . . . . . . 8.5.1 Laboratorio . . . . . . . . . . . . . . 9 Interazioni tra atomi 9.1 Approssimazione di Born-Oppenheimer . . 9.2 Superficie di Energia Potenziale . . . . . . . 9.3 Molecole biatomiche . . . . . . . . . . . . . 9.4 Solidi cristallini . . . . . . . . . . . . . . . . 9.4.1 Condizioni al Bordo Periodiche . . . 9.4.2 Teorema di Bloch . . . . . . . . . . . 9.4.3 Il potenziale vuoto . . . . . . . . . . 9.4.4 Soluzione per ll potenziale cristallino 9.4.5 Base di onde piane . . . . . . . . . . 9.5 Programma: periodicwell . . . . . . . . . .

iii

. . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

78 78 79 80 80 80 82 82 83 84 85 86 88 90 91 91 91 92 94 94 95 95 97 98 99 100 101 102 102 103

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . .

104 . 104 . 105 . 106 . 107 . 108 . 108 . 109 . 110 . 111 . 111

9.5.1

Laboratorio . . . . . . . . . . . . . . . . . . . . . . . . . . 112

A Postulati e formalismo della meccanica quantistica A.1 Stato del sistema . . . . . . . . . . . . . . . . . . . . A.2 Osservabili e misura . . . . . . . . . . . . . . . . . . A.3 Osservabili compatibili e non . . . . . . . . . . . . . A.4 Rappresentazioni . . . . . . . . . . . . . . . . . . . . A.5 Rappresentazione di Schrödinger . . . . . . . . . . . A.6 Evoluzione temporale . . . . . . . . . . . . . . . . . A.7 Definizione generale di momento angolare . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

113 . 113 . 113 . 115 . 115 . 116 . 118 . 119

B Formule utili 121 B.1 Trasformate di Legendre . . . . . . . . . . . . . . . . . . . . . . . 121 B.2 Gaussiane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 B.3 Esponenziali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 C Algoritmi utili C.1 Ricerca degli zeri . . . . . . . . . . C.1.1 Metodo di bisezione . . . . C.1.2 Metodo di Newton-Raphson C.1.3 Metodo della secante . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

123 123 123 124 124

D Software utile 125 D.1 Compilatori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 D.2 Gnuplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 D.3 BLAS e LAPACK . . . . . . . . . . . . . . . . . . . . . . . . . . 126

iv

Prefazione Lo scopo del corso di Meccanica Quantistica nell’ambito del Corso di Laurea in Fisica Computazionale è quello di trasmettere gli elementi di base necessari alla risoluzione numerica di semplici problemi di meccanica quantistica (non relativistica). E’ in pratica una sorta di laboratorio computazionale della fisica della materia, in cui vengono delineati gli ingredienti di base utilizzati nei grossi calcoli delle propriet` a di materiali, sia statici che dinamici (ossia con nuclei in movimento), che costituiscono una delle branche pi` u importanti della fisica computazionale di oggi, e con notevole rilevanza per applicazioni di interesse tecnologico e industriale. I programmi forniti nell’ambito del corso sono degli spunti. Ci si aspetta che lo studente li analizzi, li faccia girare sotto varie condizioni, studi il loro comportamento al variare dell’input, e soprattutto cerchi sempre di interpretare il loro output dal punto di vista fisico. Molti homeworks chiederanno di modificare questi programmi, aggiungendo o modificando determinate funzionalità. Il software fornito è scritto in Fortran 90. Questo linguaggio sofisticato e complesso offre gestione dinamica della memoria, operazioni su array (vettori e matrici), possibilit` a di modularizzazione dei codici, pur mantenendo una larga compatibilit` a con il Fortran 77 e senza comprometterne l’efficienza. E’ bene ricordare che la fisica computazionale è nata ben prima che fossero inventati C, Pascal e altri linguaggi, per cui esiste una enorme mole di codici e librerie scritte in Fortran 77. Dato che il Fortran 90 non è molto noto, non ci sono obiezioni all’uso da parte dello studente di un altro linguaggio come il C se vi si sente pi` u a suo agio. Una versione C di tutti i programmi è disponibile. Per eventuali approfondimenti sulla teoria della Meccanica Quantistica, esistono molti buoni libri di testo: per esempio il Griffith, lo Schiff, o i grandi classici Landau e Dirac. Per approfondimenti sull’aspetto pi` u prettamente computazionale, si rimanda a testi quali il Thijssen, il Vesely, il Koonin.

1

Bibliografia D. J. Griffiths, Introduction to Quantum Mechanics, Cummings, 2004. L. I. Schiff, Quantum mechanics, McGraw-Hill, 1968. L. D. Landau e L. Lifˇsic, Fisica Teorica, vol.3: Meccanica Quantistica non Relativistica, Editori Riuniti, 1974. J. M. Thijssen, Computational Physics, Cambridge University Press, Cambridge, 1999. Vedere anche la pagina web dell’autore: http://www.tn.tudelft.nl/tn/People/Staff/Thijssen/comphy.html. F. J. Vesely, Computational Physics - An Introduction: Second Edition, Kluwer, 2001. Vedere anche le pagine dell’autore: http://www.ap.univie.ac.at/users/Franz.Vesely/cp0102/serious.html, contenenti parti del materiale del libro. S. E. Koonin e D. C. Meredith, Computational physics - Fortran Version, Addison-Wesley, 1990. Vedere anche la pagina web di Dawn Meredith: http://pubpages.unh.edu/˜dawnm/.

2

Chapter 1

Meccanica classica Questo capitolo contiene un breve richiamo della formulazione Lagrangiana e Hamiltoniana della meccanica, seguita da una semplice applicazione numerica della meccanica classica. La prima parte ha lo scopo di richiamare dei concetti di meccanica analitica che risultano assai utili in meccanica quantistica. La seconda parte vuole mostrare come si risolvono in pratica le equazioni del moto della meccanica classica con un algoritmo di integrazione numerica. Per approfondire l’argomento, in particolare la prima parte, si consiglia: H. Goldstein, Meccanica Classica, Zanichelli, 1980 L. D. Landau e L. Lifˇsic, Fisica Teorica, vol.1: Meccanica, Editori Riuniti, 1974.

1.1

Formulazione Lagrangiana e Hamiltoniana della meccanica

La meccanica classica pu` o essere riformulata in modi alternativi, perfettamente equivalenti alla formulazione ”tradizionale” basata sull’equazione di Newton, ma pi` u comodi, eleganti e potenti. In particolare, la formulazione Hamiltoniana della meccanica introduce metodi e concetti che hanno un corrispettivo ed un’estensione naturale in meccanica quantistica. Il principio di Hamilton afferma che il moto del sistema, fra gli istanti t1 e t2 è tale per cui l’integrale di linea Z t2

I=

Ldt

(1.1)

t1

assume un valore estremo (minimo o massimo) in corrispondenza della traiettoria del moto. La funzione L ≡ L(qi , q˙i ) è detta Lagrangiana ed è funzione delle coordinate generalizzate qi e delle rispettive derivate rispetto al tempo q˙i . Le coordinate generalizzate descrivono il sistema e sono legate alle normali coordinate spaziali da una legge di trasformazione: ri ≡ ri (q1 , .., qn , t),

r˙ i ≡

X ∂ri j

∂qj

q˙j +

∂ri . ∂t

(1.2)

La Lagrangiana pu` o in generale dipendere esplicitamente dal tempo ma nel seguito assumiamo che non lo faccia.

3

Per un sistema conservativo (in cui cioè le forze derivano da un potenziale) vale L = T − V , dove T e V sono l’energia cinetica e potenziale rispettivamente. Le qi e q˙i insieme determinano lo spazio delle fasi, ovvero tutti i possibili stati del sistema. Le coordinate generalizzate in molti casi coincidono con le coordinate usuali, e le loro derivate con le velocità usuali. Dal principio di Hamilton derivano le equazioni di Lagrange (e viceversa): d dt

∂L ∂ q˙i

−

∂L = 0, ∂qi

(1.3)

che determinano il moto del sistema. Ovviamente queste sono equivalenti alle equazioni di Newton. Nella formulazione Hamiltoniana si preferisce esprimere il moto del sistema in termini delle qi e dei corrispondenti momenti generalizzati, pi , definiti come ∂L . ∂ q˙i

pi =

(1.4)

Tramite una trasformazione matematica nota come trasformazione di Legendre (vedi appendice) si introduce l’Hamiltoniana del sistema H ≡ H(pi , qi ): H(pi , qi ) =

X

q˙i pi − L(qi , q˙i ),

(1.5)

i

che determina il moto tramite le equazioni di Hamilton: p˙i = −

∂H(pi , qi ) . ∂qi

q˙i =

∂H(pi , qi ) ∂pi

(1.6)

L’Hamiltoniana altri non è che l’energia del sistema: H = T +V . In particolare, per un sistema unidimensionale, l’Hamiltoniana è H(p, x) =

p2 + V (x) ≡ T + V, 2m

(1.7)

il momento p altri non è che la quantit` a di moto (spesso chiamata anche impulso): p = mq. ˙ E’ facile verificare che le equazioni di Hamilton coincidono con quelle di Lagrange e di Newton. Da notare come coordinate e momenti (collettivamente indicate come variabili canoniche) siano considerate come variabili indipendenti, legate fra di loro solo tramite le equazioni di Hamilton. Meno usata della formulazione Lagrangiana in meccanica classica, la formulazione Hamiltoniana ha preso la sua rivincita in meccanica quantistica. In quest’ultima, l’Hamiltoniana da funzione che era diventa un onnipresente operatore e di conseguenza cambia sesso e perde la maiuscola, diventando semplicemente ”l’hamiltoniano”. Inoltre, il formalismo hamiltoniano introduce in meccanica classica una quantit` a, le parentesi di Poisson, il cui corrispettivo quanto-meccanico (il commutatore) è di importanza fondamentale. Le parentesi di Poisson [f, g] fra due funzioni f (qi , pi ) e g(qi , pi ) delle variabili canoniche sono definite come X ∂f ∂g ∂g ∂f − (1.8) [f, g] ≡ ∂qi ∂pi ∂qi ∂pi i

4

(la notazione è ovviamente mirata a mettere in evidenza la parentela con i commutatori) e godono delle seguenti proprietà: [f, f ] = 0,

[g, f ] = −[f, g]

[f, c] = 0,

(1.9)

(dove c è un numero che non dipende dalle qi e pi ) [qi , qj ] = [pi , pj ] = 0, [f, qi ] = −

∂f , ∂pi

[qi , pj ] = δij

[f, pi ] =

∂f . ∂qi

(1.10) (1.11)

Le equazioni di Hamilton possono essere espresse tramite le parentesi di Poisson: q˙i = [qi , H] ,

p˙i = − [pi , H] ,

(1.12)

come pure la derivata temporale di una funzione delle variabili canoniche: df = [f, H] dt

(1.13)

Notare che questa esprime la dipendenza temporale tramite le equazioni del moto, da non confondersi con la dipendenza esplicita dal tempo (assente in questo caso nella funzione f ).

1.2

Un’applicazione classica: moto di un punto materiale

Prima di iniziare a discutere applicazioni di meccanica quantistica, vogliamo presentare un esempio di risoluzione per via numerica di un semplice problema di meccanica classica. Lo scopo di questo esercizio è iniziare a familiarizzarci fin da subito con le tecniche per tradurre equazioni differenziali in codici di calcolo, appoggiandoci su di un esempio la cui teoria senz’altro conosciamo bene. Vogliamo ottenere la legge del moto di un punto materiale di massa m in una dimensione, soggetto ad un potenziale V (x). L’equazione differenziale che governa il moto del punto (ossia fornisce x(t)) date la posizione e la velocit` a iniziali, è la seconda legge di Newton −

d2 x dV =m 2 dx dt

(1.14)

La soluzione di questa equazione è facile da ottenere analiticamente per speciali forme di V (x), come ad esempio nel caso dell’oscillatore armonico; ma in generale ottenere una soluzione analitica potrebbe risultare assai laborioso, o impossibile. Ad esempio, V (x) stesso potrebbe non essere dato in forma analitica, ma solo in forma di una tabella numerica. Ma, soprattutto, una volta che siamo in grado di sviluppare un metodo numerico, potremo estenderlo senza troppa difficoltà a casi pi` u complessi e di interesse pratico molto maggiore, come ad esempio un sistema di molti punti materiali interagenti in tre dimensioni attraverso interazioni di coppia, o interazioni pi` u complicate.

5

1.2.1

L’algoritmo di Størmer-Verlet

La strategia generale per integrare la (1.14) è quella di suddividere l’intervallo temporale di interesse [0, T ] in N intervallini di ampiezza ∆t, sufficientemente piccoli da non commettere grossi errori approssimando la soluzione x(t) con il suo sviluppo in serie di Taylor fino ad un ordine relativamente basso, e integrare una equazione alle differenze finite per ottenere xn = x(tn ), dove tn = n∆t, n = 0 . . . N. Sviluppando in serie di Taylor attorno a xn in entrambe le direzioni: xn−1 = xn − x˙ n ∆t + (1/2)¨ xn (∆t)2 − (1/6)x ¨˙ n (∆t)3 + O[(∆t)4 ] 2 xn+1 = xn + x˙ n ∆t + (1/2)¨ xn (∆t) + (1/6)x ¨˙ n (∆t)3 + O[(∆t)4 ]

(1.15)

(la notazione x˙ indica la derivazione rispetto al tempo) e sommandole fra loro si ottiene (1.16) xn+1 = 2xn − xn−1 + x ¨n (∆t)2 + O[(∆t)4 ] Ora utilizziamo la legge di Newton (1.14), o

1 dV x ¨n = − ≡ fn m dx x=xn

(1.17)

xn+1 = 2xn − xn−1 + fn (∆t)2 + O[(∆t)4 ].

(1.18)

ottenendo Questa equazione ci permette, nota la posizione ai tempi n − 1 e n e la forza al tempo n, di ottenere una stima della posizione al tempo n + 1, e quindi fornisce un algoritmo (detto algoritmo di Størmer-Verlet) per ottenere iterativamente la traiettoria x(t) del punto sotto forma di una tabella numerica.

1.3

Programma: newton

Il programma newton.f901 (oppure la sua versione C, newton.c2 ) implementa l’algoritmo di Størmer-Verlet per un punto materiale soggetto ad un potenziale V (x). Nel programma è definito in particolare 2π (1.19) 5 il cui andamento nella regione |x| ≤ 6 è indicato in fig. 1.1. Si tratta di un potenziale parabolico (come se fosse quello di un oscillatore armonico), ma ”modulato” attraverso un termine oscillante con una periodicità ` solo un esempio, pari a 5 unit` a di lunghezza. Non è un potenziale famoso! E scelto in modo assolutamente arbitrario, di un potenziale con pi` u posizioni di equilibrio, per rendere le simulazioni pi` u interessanti. Il potenziale è definito in una subroutine separata nel programma, e può essere variato facilmente senza dover modificare il corpo principale contenente l’algoritmo di integrazione. Per semplicit` a si assume inoltre una massa unitaria. E’ facile vedere che definire una massa diversa è equivalente a moltiplicare il time step ∆t per un √ fattore m, e quindi tale assunzione non comporta alcuna perdita di generalità. V (x) = x2 [2 − cos(kx)]

1 2

,

k=

http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/F90/newton.f90 http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/C/newton.c

6

70

V(x)

60

50

40

30

20

10

0

-6

-4

-2

0

2

4

6

Figure 1.1: Potenziale V(x) definito in newton.

1.3.1

Struttura del programma

Come indicato sopra, il potenziale V (x) è definito dalla subroutine potential separata, e quindi è in un certo senso esterno al programma principale. Per un dato x, potential ritorna il potenziale V (x) e la corrispondente forza F (x) = −dV /dx. Solo quest’ultima è, strettamente parlando, necessaria per portare avanti l’integrazione; ma conoscere anche il potenziale permette di verificare la legge di conservazione dell’energia passo per passo. Il programma chiede i seguenti dati di input: time step ∆t, posizione iniziale x1 , energia totale E1 , numero di time steps (verrà seguita l’evoluzione temporale per un tempo totale N ∆t), e nome del file di uscita (è anche possibile convogliare l’output sullo standard output). L’energia totale serve a calcolare la velocità iniziale v1 , secondo 1 E1 = mv12 + V (x1 ) 2

(1.20)

Naturalmente è necessario che sia E1 ≥ V (x1 ) (se cos`ı non è il programma si ` stato scelto di usare l’energia, piuttosto che direttamente la velocità arresta). E iniziale, perchè cos`ı è immediatamente chiaro quale regione del potenziale verrà esplorata nel corso della simulazione: quella in cui V (x) ≤ E1 . I punti in cui V = E1 sono i punti di inversione, in cui la velocità è nulla e la particella inverte la sua direzione di moto. Inoltre, come ben noto l’energia totale 1 E = mv 2 + V (x) 2

(1.21)

è supposta essere una costante del moto. La conservazione dell’energia costituisce per questi problemi un eccellente strumento di verifica che non vi siano errori grossolani o siano stati dati input errati, come un ∆t troppo grande. Assegnare la velocit` a iniziale significa, nell’ottica dell’algoritmo di StørmerVerlet, assegnare la posizione al ”tempo precedente”. Tale algoritmo però opera su posizioni ed accelerazioni. La velocità alla posizione n può essere in generale

7

ottenuta come vn =

xn+1 − xn−1 2∆t

(1.22)

con un errore dell’ordine di (∆t)2 (come subito si vede dalle (1.15), sottraendo un’equazione dall’altra). Tuttavia, questa formula coinvolge tre punti e come condizioni iniziali noi abbiamo a disposizione la sola posizione e velocità iniziali al medesimo punto x1 . Si è allora scelto di utilizzare l’espressione v1 '

x1 − x0 ∆t

(1.23)

per definire la ”posizione precedente” x0 a partire da x1 e v1 . Se x0 e x1 fossero dati, ad una tale formula sarebbe associato un errore dell’ordine di ∆t nella velocit` a, quindi non si tratta di una buona stima. Tuttavia, quello che al pi` u succede è che una volta che la simulazione è partita, l’energia totale si assesterà ad un valore lievemente diverso da quello da noi richiesto. In uno schema pi` u raffinato, si potrebbe pensare di effettuare una correzione al secondo passo di integrazione per ottenere una migliore corrispondenza fra energia richiesta e energia effettiva. La velocit` a iniziale è sempre definita positiva dal programma, ossia il punto inizia a muoversi verso gli x positivi. Data la simmetria del potenziale, questa non è una limitazione (invertire la velocità è equivalente a lasciarla invariata e cambiare il segno della coordinata di posizione iniziale). Il loop sul numero di time steps è la parte centrale del programma. Il punto su cui va posta particolare attenzione è il trasferimento dell’informazione della ”posizione precedente” da una iterazione alla successiva, effettuato usando tre variabili x, xprev, xprevsav. A fini di controllo sulla conservazione dell’energia, viene calcolata anche la velocit` a e l’energia cinetica. Per ottenere una precisione accettabile, le velocità vengono definite utilizzando la (1.22). Questo però comporta che per definire vn è necessario conoscere xn+1 , e quindi che alla conclusione di ogni iterazione le velocit` a calcolate sono riferite all’iterazione precedente, ossia alla posizione che aveva la particella all’inizio dell’iterazione, prima di essere spostata. I dati emessi sul file di output sono, nell’ordine, il numero dell’iterazione, la posizione, la velocit` a, la forza, l’energia cinetica, l’energia potenziale e l’energia totale: tutte riferite alla posizione che aveva la particella all’inizio dell’iterazione.

1.3.2

Laboratorio

Ecco alcuni spunti per esperimenti numerici che si possono fare con questo codice. Un valore ”sicuro” per il time step è attorno a 0.01. • Graficare le varie quantit` a in funzione del tempo ed interpretarle. • Graficare le quantit` a tra loro (molto semplice con gnuplot) (vedere Appendice D.2), ad esempio la velocità in funzione della posizione (è una vista nello spazio delle fasi).

8

• Confrontare il risultato di un input 0.005/0/28.36/3000 con un input 0.005/0/28.35/3000 (un’ispezione del potenziale può anticipare il risultato...) • Provare ad aumentare il ∆t. Tenere sotto controllo la conservazione dell’energia e le traiettorie nello spazio delle fasi. Discutere i risultati con l’input 0.2/0/27/10000. Ed ecco inoltre spunti per modifiche del codice: • Definire altri potenziali. • Immagazzinare il potenziale in una tabella numerica, e effettuare interpolazione in questa tabella durante l’integrazione (molto utile nel caso in cui il calcolo del potenziale sia dispendioso: in questo modo basta farlo una volta per tutte). • Introdurre una forza di attrito (F = −γ x). ˙ • Introdurre una seconda massa identica alla prima, ed assumere che le due masse siano accoppiate tra loro, ad esempio in modo armonico: V (x2 − x1 ) = (1/2)k(x2 − x1 )2 , e che possano ”incrociarsi”. • Introdurre una seconda massa identica alla prima, e lasciare che si muovano indipendentemente, facendole tuttavia collidere elasticamente quando entrano in contatto.

9

Chapter 2

Introduzione all’equazione di Schr¨ odinger Questo capitolo ci dar` a una prima presa di contatto con la meccanica quantistica ”vera”, e in particolare introdurrà l’oggetto matematico che il pi` u delle volte finiremo con il calcolare per ottenere le proprietà di un sistema microscopico: la funzione d’onda, soluzione dell’equazione di Schrödinger.

2.1

Verso la Meccanica Quantistica

Ricapitoliamo qui, senza nessuna pretesa di accuratezza storica, i principali risultati che portarono fra il 1900 e il 1930 all’introduzione della Meccanica Quantistica nella forma che conosciamo adesso.

2.1.1

Dualismo onda-particella: fotoni

L’effetto fotoelettrico, lo spettro di radiazione del corpo nero e l’effetto Compton dimostrano che la radiazione elettromagnetica di frequenza ν si comporta come se fosse costituita da fotoni, particelle di massa nulla, energia E = hν e quantità di moto p = E/c; la costante h ha le dimensioni di un’azione (energia per tempo, o quantit` a di moto per posizione) ed è detta costante di Planck. Tuttavia queste particelle continuano a comportarsi sotto certi aspetti come onde, da cui il cosidetto dualismo onda-particella. Entrambe gli aspetti, corpuscolare e ondulatorio, sembrano presenti allo stesso tempo ed entrambe necessari per spiegare il comportamento della radiazione. Esperimenti di interferenza con fenditure sono particolarmente istruttivi e sorprendenti a questo riguardo. In un tipico setup sperimentale, una sorgente luminosa monocromatica emette un fascio di luce che colpisce uno schermo con due fenditure; la luce attraversa le fenditure e colpisce un altro schermo piazzato ad una certa distanza dietro al primo. Su questo si osservano le cosidette frange di interferenza. La spiegazione del fenomeno in termini di onde è banale: a seconda della differenza di cammino ottico, le onde provenienti dalle due fenditure si sommano in fase o in controfase, dando origine a zone chiare e scure alternate. La spiegazione in termini di fotoni è quantomeno problematica,

10

tanto pi` u che si possono mostrare (sperimentalmente!) due aspetti notevoli del problema: • l’interferenza si realizza anche se l’intensità della luce è cos`ı debole da avere (in media) un solo fotone presente sul percorso fra la sorgente e lo schermo su cui si osservano le frange di interferenza; • l’interferenza sparisce se si tappa una delle due fenditure, o anche se si cerca in qualunque modo di misurare ”da quali delle due fenditure è passato il fotone”. E’ quindi necessario assumere che un fotone ”interferisca con se stesso”, ma solo se non lo costringiamo a rivelare la sua natura corpuscolare con una ”misura” che ci dica ”da dove è passato”. Ovvero: natura ondulatoria e corpuscolare del fotone sembrano inestricabilmente legate.

2.1.2

Quantizzazione dei livelli di energia

Il problema della stabilit` a dell’atomo e della forma degli spettri atomici ricevettero una prima soluzione con la cosiddetta ”vecchia meccanica quantistica” di Niels Bohr. Bohr postul` o che l’atomo fosse stabile solo per certi valori discreti dell’energia; l’assorbimento o l’emissione di un fotone di frequenza ν potevano avvenire solo in seguito ad una transizione fra tali valori discreti (livelli di energia) in cui la differenza di energia fosse pari ad hν. Nel modello di Bohr, le orbite permesse nell’atomo sono quelle per cui vale la relazione Z

l=

pdq = nh

(2.1)

orbita

ovvero il momento angolare ha solo valori discreti, multiplo della costante di Planck. Tale idea pu` o essere estesa ad altri problemi semplici: l’oscillatore armonico, la buca di potenziale, etc. Il concetto di quantizzazione dei livelli di energia e la presenza di livelli discreti risolve inoltre anche molti enigmi sui calori specifici (un altro vecchio problema della meccanica classica). La vecchia meccanica quantistica si dimostrò quindi un’intuizione vincente, ma rimaneva una teoria basata su assunzioni ad hoc e di difficile estensione a sistemi un minimo complicati.

2.1.3

Dualismo onda-particella: elettroni

La presenza della costante di Planck sia nella relazione energia-frequenza dei fotoni che nella condizione di Bohr suggerisce che h sia una costante universale e fondamentale non solo per i fotoni ma anche per le particelle. Inoltre la condizione di Bohr pu` o essere riscritta sotto una forma che è molto suggestiva della presenza di un carattere ondulatorio anche nelle particelle: per orbite circolari di circonferenza L, pL = nh, ovvero L = nh/p = nλ. La lunghezza λ ha la stessa definizione della lunghezza d’onda di un fotone di quantità di moto p. La condizione di Bohr diventa quindi come la condizione di risonanza per le onde: il numero di ”lunghezze d’onda” che stanno su di una circonferenza è

11

intero. De Broglie sugger`ı che in generale si può associare un’onda di lunghezza d’onda λ = h/p ad una particella di quantità di moto p. Il significato preciso di tale ”associazione” sar` a chiarito in seguito. Il dualismo onda-particella prende cos`ı un aspetto ”simmetrico”: non solo le onde hanno un carattere di particella, ma anche le particelle hanno un carattere di onda. La prova finale del carattere ondulatorio associato alle particelle è data dall’osservazione di un fenomeno tipicamente ondulatorio quale la diffrazione di Bragg in fasci di particelle. Questo fu realizzato nel 1927 da Davisson e Gerner con fasci di elettroni. Adesso sappiamo che la diffrazione di Bragg di particelle è osservabile persino con fasci di atomi di Elio.

2.1.4

Principio di indeterminazione

Una conseguenza del dualismo onda-particella è una limitazione intrinseca alla precisione delle misure, in particolare alla possibilità di misurare contemporaneamente quantit` a di moto e posizione di una particella con precisione arbitraria. Per determinare la seconda serve osservare la particella con luce di lunghezza d’onda sempre pi` u piccola; ma questa avrà un’energia associata sempre pi` u grande e perturber` a la quantità di moto della particella. Visto in un altro modo: se le particelle hanno anche un carattere ondulatorio, le si possono pensare come dei ”pacchetti d’onda”, sovrapposizione di onde monocromatiche; ma un pacchetto d’onda tanto pi` u è localizzato nello spazio, tanto pi` u ampio è nello spazio reciproco (che è poi quello delle quantità di moto, come vedremo). Heisenberg formalizz` o questo aspetto nel principio, o pi`‘u esattamente, relazione di indeterminazione: per il prodotto fra l’incertezza, sulla posizione ∆x e sulla quantit` a di moto ∆px , vale la diseguaglianza ∆x∆px ≥ ¯h/2, dove h ¯ = h/2π.

2.2

L’equazione di Schr¨ odinger per una particella

Abbiamo visto nella sezione precedente quale fenomenologia la Meccanica Quantistica debba poter spiegare. Un passo avanti decisivo fu compiuto da Schrödinger, che estese e formalizz` o l’idea di ”onda associata ad una particella” con l’introduzione della funzione d’onda. La funzione d’onda è determinata come soluzione di un’equazione differenziale, l’equazione di Schr¨ odinger, che ammette soluzione solo per determinati valori dell’energia, corrispondenti ai livelli energetici della particella. Il quadrato della funzione d’onda determina la probabilità di trovare la particella in una determinata posizione, in seguito ad una misura. La posizione della particella è quindi descritta in modo probabilistico: data una funzione d’onda, sappiamo qual è la probabilità che la particella sia in un certo punto, ma non ”dove sta esattamente”. Consideriamo per prima l’equazione di Schrödinger indipendente dal tempo per un sistema particolarmente semplice: una particella in un potenziale unidimensionale V (x). L’equazione ha la forma seguente: −

¯h2 d2 ψ + V (x)ψ(x) = Eψ(x), 2m dx2

12

(2.2)

dove ψ(x) è la funzione d’onda (in generale complessa). Da dove scappa fuori? Uno dei postulati della meccanica quantistica è che la quantità di moto lungo la direzione x è rappresentata dall’operatore differenziale: p = −i¯h

d dx

(2.3)

Possiamo quindi riscrivere la (2.2) come "

#

p2 + V (x) ψ(x) = Eψ(x). 2m

(2.4)

L’operatore fra parentesi quadre a sinistra rappresenta quindi una vecchia conoscenza: è l’Hamiltoniano del sistema, H = p2 /2m + V . Ciò ci permette di identificare E come l’energia totale (cinetica+potenziale) del sistema. L’equazione non ha necessariamente soluzioni fisicamente accettabili (ovvero non divergenti) per qualunque valore di E: i valori di di E per i quali una soluzione esiste rappresentano quindi i livelli di energia ammessi. Notiamo come nella (2.4) l’operatore H applicato alla soluzione ψ(x) produce la ψ(x) medesima, moltiplicata per una costante. Questo si riflette nella seguente terminologia: la soluzione dell’Eq.(2.2) è detta autofunzione dell’Hamiltoniano, il corrispondente valore di E è detto autovalore. Nei prossimi capitoli vedremo in modo pi` u chiaro il perché di questa terminologia proveniente dall’algebra lineare. La funzione d’onda è definita a meno di una costante arbitraria. Si usano di solito funzioni d’onda normalizzate: Z

|ψ(x)|2 dx = 1.

(2.5)

L’energia E è il valore di aspettazione dell’operatore H sullo stato definito da ψ, cioè: Z

E=

"

#

p2 ψ (x) + V (x) ψ(x)dx = 2m ∗

Z

ψ ∗ (x)Hψ(x)dx.

(2.6)

Per funzioni d’onda non normalizzate: R

E=

ψ ∗ (x)Hψ(x)dx R . |ψ(x)|2 dx

(2.7)

In generale, quanto sopra si applica a qualunque operatore f (x, p), in cui p è sostituito dall’operatore differenziale (refqdm). Si noti che H è un operatore: un oggetto che agisce su di una funzione e ha come risultato un’altra funzione, in questo caso. Non è lecito spostare la ψ da destra a sinistra di H o viceversa !!

2.2.1

L’equazione di Schr¨ odinger dipendente dal tempo

In generale, l’equazione di Schr¨ odinger deve tenere conto anche della dipendenza temporale. L’equazione dipendente dal tempo ha la forma: −

¯h2 ∂ 2 Ψ(x, t) ∂Ψ(x, t) + V (x, t)Ψ(x, t) = i¯h 2 2m ∂x ∂t

13

(2.8)

dove Ψ(x, t) è la funzione d’onda. Nel caso importante in cui il potenziale non dipenda esplicitamente dal tempo: V (x, t) = V (x), ci si può ricondurre alle soluzioni dell’equazione indipendente dal tempo, Eq.(2.2): vediamo come. Scriviamo la Ψ(x, t) come prodotto di una funzione di solo x e una di solo t: Ψ(x, t) = ψ(x)f (t)

(2.9)

Sostituendo la (2.9) nella (2.8) e dividendo per ψ(x)f (t) si trova −

1 ¯h2 d2 ψ(x) i¯h df + V (x) = 2 ψ(x) 2m dx f (t) dt

(2.10)

In questa equazione il membro sinistro dipende solo da x, e quello destro solo da t. Entrambi i membri devono allora essere uguali ad una costante, che chiamiamo E. Otteniamo cos`ı due equazioni: quella spaziale è appunto la (2.2), mentre quella temporale è iE df = − f (t) dt ¯h

(2.11)

f (t) = Ce−iEt/¯h

(2.12)

la cui soluzione è banalmente

dove C è una costante (determinata dalla normalizzazione). Si tratta di un punto nello spazio complesso in rotazione attorno all’origine con frequenza angolare E/¯h. Abbiamo quindi trovato che ad ogni soluzione ψn (x) della (2.2), corrispondente a un certo valore di En , corrisponde anche una soluzione della (2.8) Ψn (x, t) = ψn (x)e−iEn t/¯h

(2.13)

Questa soluzione è uno stato stazionario, perchè |Ψn (x, t)|2 = |ψn (x)|2 non dipende dal tempo. Ψn (x, t) data dalla (2.13) non è una soluzione generale dell’equazione di Schrödinger dipendente dal tempo (2.8). Tuttavia, si può dimostrare che l’insieme di tutte le soluzioni possibili è costituito dalle combinazioni lineari delle autofunzioni dell’energia, ossia qualsiasi soluzione Ψ(x, t) è sempre esprimibile come una sovrapposizione di stati stazionari: Ψ(x, t) =

X

cn Ψn (x, t).

(2.14)

n

Questo è un risultato importante: data una funzione d’onda Ψ(x, t0 ) che si sa essere una soluzione valida ad un certo istante t0 , la sua evoluzione temporale può essere ottenuta facilmente se si riesce a svilupparla in stati stazionari al tempo t0 secondo la (2.14). La soluzione numerica diretta [ossia operando nello spazio (x, t)] dell’equazione (2.8) è in generale un problema difficile, che porta spesso a instabilità numeriche. Quasi sempre l’evoluzione temporale di una funzione d’onda non corrispondente ad uno stato stazionario viene perciò studiata decomponendola in autofunzioni dell’energia—la cui evoluzione temporale è data dalla (2.13)—secondo la (2.14). Tuttavia, nel caso generale in cui il potenziale dipende dal tempo la separazione delle variabili non è possibile, e il problema va quindi affrontato risolvendo direttamente la (2.8).

14

2.3

La particella libera

Consideriamo una particella libera in una dimensione, la cui equazione di Schrödinger dipendente dal tempo è (dalla (2.8) per V (x) = 0): ∂Ψ(x, t) ¯h2 ∂ 2 Ψ(x, t) = i¯h (2.15) 2m ∂x2 ∂t La soluzione di questa equazione, come facilmente si vede, sono le onde piane: −

Ψ(x, t) = Cei(±kx−ωt)

(2.16)

dove C è una costante (scelta in modo da normalizzare la funzione correttamente), e k e ω sono fra loro legati dalla relazione ¯h2 k 2 =h ¯ω (2.17) 2m E può assumere qualsiasi valore reale positivo. Non vi è quantizzazione: l’energia di una particella libera pu` o avere qualsiasi valore. Notare come esistano due diverse funzioni d’onda per ogni valore di E: nel gergo della Meccanica Quantistica, si dice che ogni autovalore è due volte degenere. Come sappiamo, potevamo anche considerare l’equazione di Schrödinger indipendente dal tempo per lo stesso problema E=

−

¯h2 d2 ψ = Eψ(x) 2m dx2

(2.18)

ψ(x) = Ce±ikx

(2.19)

e ottenere una soluzione dove E =

h ¯ 2 k 2 /2m,

e poi dire, applicando la (2.13), che Ψ(x, t) = ψ(x)e−iEt/¯h

(2.20)

Il risultato è chiaramente identico. La (2.16) e la (2.20) rappresentano un’onda che si propaga con velocità v = ω/k = h ¯ k/2m. La loro forma può lasciare perplessi. In primo luogo, la funzione è delocalizzata in modo uniforme su tutto lo spazio: |Ψ(x, t)| = 1. Osserviamo per` o che la (2.16), di cui scegliamo la soluzione con segno positvo. è un’autofunzione dell’operatore quantità di moto, definito in (2.3), con autovalore ¯hk: d (2.21) pΨ(x, t) ≡ −i¯h Ψ(x, t) = ¯hkΨ(x, t). dx Per analogia con l’equazione di Schrödinger, ciò implica che con questa funzione d’onda la quantit` a di moto ha un valore ben definito ¯hk, e che quindi la sua indeterminazione è nulla: ∆p = 0. Torneremo nel Cap.A su questo aspetto. La relazione di indeterminazione ci dice allora che la posizione è completamente indeterminata: ∆x = ∞. Non è quindi preoccupante che la nostra funzione d’onda Rsia uniformemente delocalizzata; tuttavia ciò la rende non normalizzabile: |Ψ(x, t)|dx → ∞. Dovremo abituarci a convivere con oggetti la cui corretta definizione matematica passa per una procedura di limite. In questo caso, potremmo considerare la particella libera come limite per L → ∞ di una particella in una scatola di dimensioni L, con la condizione ψ(x + L) = ψ(x) (condizioni periodiche ai bordi).

15

2.4

Pacchetti d’onda

Come abbiamo visto, la soluzione (2.16) per la particella libera non somiglia molto al moto di una particella libera classica, in quanto: 1. l’ampiezza della (2.16) è costante; 2. se la quantit` a di moto classica p deve corrispondere a h ¯ k, la velocità dell’onda sembra essere la metà di ciò che ci si aspetterebbe. Per poter ottenere un limite classico sensato, dobbiamo introdurre il concetto di pacchetto d’onde, e assumere che il moto classico si ottenga sommando fra loro molte onde piane del tipo (2.16), anzichè considerando un’onda sola. Proviamo a considerare per un momento la sola parte spaziale (non è una limitazione: è la soluzione della corrispondente equazione di Schrödinger indipendente dal tempo, e conosciamo il suo legame (2.13) con la soluzione completa), e ipotizziamo una soluzione oscillante come un’onda piana con un certo vettore d’onda k0 , ma localizzata nello spazio in una regione di lunghezza L: ψ(x) = eik0 x se |x| ≤ L/2 = 0 se |x| > L/2

(2.22)

Ci chiediamo se la (2.22) è una soluzione dell’equazione (2.18). Per fare questo, appoggiamoci sulla teoria delle trasformate di Fourier, secondo cui qualunque ψ(x) pu` o essere espressa in termini di uno sviluppo in onde: Z +∞

ψ(x) =

F (k)eikx dk

(2.23)

−∞

dove le ampiezze F (k) si possono ottenere da una ψ(x) mediante una trasformata inversa, Z 1 +∞ F (k) = ψ(x)e−ikx dx (2.24) 2π −∞ Nel nostro caso, la forma particolare (2.22) che abbiamo ipotizzato dà F (k) =

1 2π

Z L/2

e−i(k−k0 )x dx =

−L/2

L sin[(k − k0 )L/2] 2π (k − k0 )L/2

(2.25)

Come noto, la funzione sin y/y ha un picco di ampiezza 1 a y = 0, si annulla per y = ±π, e presenta altre oscillazioni di ampiezza molto inferiore al picco principale, che decadono come 1/y al crescere di y. Pertanto F (k) ha un picco di altezza massima L/2π e larghezza a metà altezza approssimativamente ∆k ∼ 2π/L. Abbiamo quindi trovato che è possibile costruire una soluzione localizzata e oscillante con numero d’onda k0 , ma per fare questo dobbiamo sovrapporre un insieme di onde piane con numero d’onda centrato attorno a k0 ma con una dispersione ∆k. La (2.22) descrive pertanto una particella quantistica la cui posizione è determinata con una incertezza ∆x ∼ L, e la cui quantità di moto è determinata con una incertezza ∆p = ¯h∆k ∼ h/L. Abbiamo quindi ∆x∆p ∼ h, che è l’espressione del principio di indeterminazione. Il problema di una singola onda

16

piana è che la sua quantit` a di moto è determinata esattamente, e questo rende la posizione totalmente indefinita. Analoghe considerazioni possono essere effettuate per quanto riguarda la variabile temporale. In questo caso si effettuano trasformate di Fourier tra la variabile temporale e lo spazio delle frequenze. Se un treno d’onde ha una durata finita complessiva T (che sarà quindi il ∆t) ed effettua N oscillazioni, la precisione nella determinazione della sua frequenza è circa pari a 1 oscillazione, ossia 1 2π/ω0 ∆ω ∼ (2.26) = ω0 N T da cui ∆t∆ω ∼ 2π (2.27) ovvero l’indeterminazione nel tempo e quella nell’energia sono legate da ∆t∆E ∼ h. Il limite classico della meccanica quantistica passa quindi necessariamente attraverso i pacchetti d’onda per poter confinare la particella in una regione finita. Va notato che alla dispersione in k corrisponderà anche una dispersione in energia. Ogni componente k soddisfa all’equazione di Schrödinger indipendente dal tempo per l’energia E = ¯h2 k 2 /2m. Le componenti si sommano solo dopo aver moltiplicato ciascuna di esse per il fattore di fase dipendente dal tempo, secondo la (2.20). Questo fa s`ı che l’aspetto del pacchetto possa in generale variare nel tempo. Occupiamoci allora della seconda questione relativa alla velocità. Immaginiamo di costruire un semplice pacchetto costituito da due sole onde, una di numero d’onda k0 − δk e una di numero d’onda k0 + δk, dove δk è piccolo. Le frequenze angolari corrispondenti (attraverso la 2.17) saranno ω0 −δω e ω0 +δω. Ψ(x, t) = ei(k0 −δk)x e−i(ω0 −δω)t + ei(k0 +δk)x e−i(ω0 +δω)t = ei(k0 x−ω0 t) [2 cos(δkx − δωt)]

(2.28)

ossia un’onda piana di numero d’onda k0 modulata da un fattore oscillante con un numero d’onda assai pi` u piccolo, ossia con una lunghezza d’onda molto pi` u grande. Questo è un inviluppo analogo a quello che dà luogo ai battimenti in acustica. L’inviluppo si muove con una velocit` a diversa da quella dell’onda che contiene. Possiamo trovare la sua velocità seguendo ad esempio lo spostamento nel tempo del massimo corrispondente a un argomento nullo del coseno: δkx − δωt = 0 ossia

(2.29)

dω (2.30) dk La quantit` a vg è detta velocit` a di gruppo. Dato un pacchetto d’onde qualsiasi, per ogni coppia di componenti vicine si può pensare che valga la (2.30), che quindi rappresenta la velocit` a del pacchetto stesso. Nel limite classico, è la velocit` a di gruppo che diventa la velocità della particella classica. Dalla (2.17) si ha subito infatti ¯hk vg = (2.31) m x = vg t

,

17

vg =

che è quanto ci si aspetta. Nel caso di una particella libera, la velocità del pacchetto è quindi doppia rispetto a quella dell’onda, e corrisponde al limite classico. Da notare infine che il pacchetto si delocalizza nel tempo perchè ciascuna delle sue componenti k si propaga con una velocità ω/k diversa da quella delle altre componenti. Affinchè il pacchetto non si degradi, occorrerebbe che ω/k fosse una costante. Questo è in effetti il caso delle onde elettromagnetiche nel vuoto, ma non delle onde associate a particelle con massa finita.

2.5

Potenziali modello

Consideriamo qualche caso di potenziale semplice. I casi seguenti potranno sembrare artificiosi e di dubbio interesse, ma in realtà molti sistemi fisici sono descrivibili in modo approssimato con potenziali molto semplici.

2.5.1

Gradino di potentiale

Consideriamo il seguente potenziale a gradino: V (x) = 0 per x < 0, V (x) = W per x > 0 (W > 0). Si presenta subito una difficoltà: cosa succede alla funzione d’onda nel punto di discontinuit` a, x = 0, del potenziale? La risposta è nota dalla teoria matematica, ma possiamo darne una ”fisica” considerando il potenziale discontinuo come limite di potenziali continui che passano da V (0) = 0 a V () = W per → 0. Riscriviamo l’equazione di Schrödinger come: ψ 00 (x) =

2m (V (x) − E) ψ(x) ¯h2

(2.32)

da cui si ricava l’ovvio risultato che la derivata seconda della funzione d’onda è discontinua in x = 0. Integriamo fra x = 0 e x = : ψ 0 () − ψ 0 (0) =

2m ¯h2

Z

(V (x) − E) ψ(x)dx.

(2.33)

0

Siccome sia V (x) che ψ(x) sono finiti nell’intervallo (0, ), il secondo membro tende a 0 per → 0. Quindi ψ 0 () → ψ 0 (0). Analogamente si dimostra la continuit` a di ψ(x) intorno a x = 0. Queste sono le condizioni da imporre ovunque sia presente un gradino (finito) di potenziale. Per risolvere il problema, si devono distinguere tre intervalli di energia: 0) E < 0: non esistono soluzioni, o pi` u esattamente, le sole soluzioni sono esponenziali reali, che divergono e quindi non sono fisiche. 1) E > W : le soluzioni sono onde piane, ψ(x) = Aeikl x + Be−ikl x ,

√ x < 0,

kl =

ψ(x) = Ceikr x + De−ikr x , x > 0, kr =

2mE/¯h;

(2.34)

q

2m(E − W )/¯h.(2.35)

Le condizioni di continuit` a a x = 0 ci danno A+B = C +D

(2.36)

kr (A − B) = kl (C − D).

(2.37)

18

Abbiamo quindi quattro incognite e due condizioni, pi` u la normalizzazione (nel nostro caso, arbitraria). Rimane quindi una costante indeterminata. In effetti, ci sono due soluzioni per ogni valore di energia e quindi ogni combinazione lineare delle due è una soluzione accettabile. Possiamo per esempio selezionale la soluzione corrispondente ad un’onda incidente da sinistra e trasmessa a destra (D = 0). In questo caso, si trova √ √ B E− E−W C B √ =√ =1+ . , (2.38) A A A E+ E−W 2) 0 < E < W : le soluzioni sono onde piane per x < 0, onde evanescenti per x > 0: √ ψ(x) = Aeikl x + Be−ikl x , x < 0, kl = 2mE/¯h; (2.39) ψ(x) = Ce−kr x ,

x > 0,

kr =

q

2m(W − E)/¯h.

(2.40)

L’esponenziale con il segno opposto per x > 0 è ovviamente non accettabile! Le condizioni di continuità a x = 0 ci danno kr (A − B) = −kl C

A + B = C, da cui

√ √ B E−i W −E √ = √ , A E+i W −e

B C =1+ . A A

(2.41)

(2.42)

In questo caso c’è una sola soluzione per ogni valore di E: l’autovalore è non degenere (cosa prevista dal teorema di non degenerazione, valido in sistemi unidimensionali). Da notare come nell’intervallo 2 la funzione d’onda nella regione x > 0, classicamente inaccessibile, è evanescente e rapidamente tendente a zero, ma comunque non nulla: esiste una probabilit` a piccola ma finita di trovare la particella in una zona dove ”non dovrebbe stare”. E’ un fenomeno intrinsecamente quantomeccanico, di grande rilevanza.

2.5.2

Barriera di potenziale

Consideriamo ora una barriera di potenziale: V (x) = W per |x| < a/2, V (x) = 0 per x < −a/2 e x > a/2. In questo caso si richiedono due operazioni di ”matching” della funzione d’onda, a x = −a/2 e x = a/2. Per gli intervalli di energie 0 e 1 sopra introdotti, i risultati sono del tutto analoghi al caso del gradino: nessuna soluzione e due soluzioni degeneri per ogni E, rispettivamente. Pi` u interessante il caso dell’intervallo 2. Consideriamo una soluzione che si propaga verso destra nella regione x < −a/2. Tale soluzione diventerà un’onda evanescente nella regione classicamente proibita |x| < a/2. Tuttavia in x = a/2 sopravviver` a una componente piccola di onda evanescente che avrà come corrispettivo un’onda propagantesi nella regione x > a/2. Esistono quindi soluzioni che ”scavalcano” la barriera: è il cosiddetto effetto tunnel.

19

Scriviamo la funzione d’onda per il caso come quello descritto sopra, assumendo per semplicit` a il coefficiente dell’onda incidente uguale a 1: √ (2.43) ψ(x) = eikx + Ae−ikx , x < −a/2, k = 2mE/¯h; 0

0

ψ(x) = Be−k x + B 0 e−k x , ikx

ψ(x) = Ce

,

|x| < a/2,

k0 =

q

2m(W − E)/¯h; (2.44)

x > a/2.

(2.45)

Notare la presenza di una componente di onda crescente nella funzione d’onda per |x| < a/2: in effetti, non abbiamo il diritto di escluderlo a priori! Il calcolo, relativamente semplice ma un po’ laborioso, dà il seguente risultato: 4E(W − E) 4E(W − E) + W 2 sinh2 (k 0 a) W 2 sinh(k 0 a) = 1 − |C|2 = . 4E(W − E) + W 2 sinh2 (k 0 a)

|C|2 =

(2.46)

|A|2

(2.47)

Il caso E > W si pu` o ottenere con la sostituzione 0

p

k −→ ik1 = i

2m(E − W ) ¯h

(2.48)

e dà il seguente risultato: 4E(E − W ) 4E(E − W ) + W 2 sin2 (k1 a) W 2 sin(k1 a) = 1 − |C|2 = . 4E(E − W ) + W 2 sin2 (k1 a)

|C|2 =

(2.49)

|A|2

(2.50)

|C|2 è detto coefficiente di trasmissione, mentre |A|2 è il coefficiente di riflessione. In generale, essi dipendono dall’energia e dalla forma del potenziale.

2.5.3

Buca di potenziale

La buca di potenziale è una schematizzazione molto semplice di un potenziale attrattivo o vincolante: V (x) = −W per |x| < a/2, V (x) = 0 per x < −a/2 e x > a/2. Consideriamo l’intervallo di energia interessante: −W < E < 0. Scriviamo la soluzione sotto la forma 0

k0 =

ψ(x) = Aek x , ψ(x) = B cos(kx − α), −k0 x

ψ(x) = Ce

,

k=

q

2m|E|/¯h,

q

2m(W + E)/¯h,

x ≥ a/2,

x ≤ a/2 |x| ≤ a/2

(2.51) (2.52) (2.53)

pi` u conveniente per i calcoli (scrivere la soluzione come coseno + fase è del tutto equivalente a scrivere come somma di esponenziali complessi: possiamo sempre ricondurci a soluzioni reali). Le condizioni di continuità a x = −a/2 e x = a/2 sono: 0

Ae−k a/2 = B cos(−ka/2 − α),

0

k 0 Ae−k a/2 = −kB sin(−ka/2 − α), (2.54)

20

0

0

Ce−k a/2 = B cos(ka/2 − α),

−k 0 Ce−k a/2 = −kB sin(ka/2 − α),

(2.55)

ovvero, dividendo membro a membro, k tan(ka/2 + α) = k 0 ,

k tan(ka/2 − α) = k 0 .

(2.56)

Queste due condizioni possono essere soddisfatte contemporaneamente solo se α = 0 o se α = π/2. Distinguiamo i due casi: • Soluzioni pari (α = 0): esistono solo a energie per cui k tan(ka/2) = k 0 . Introduciamo le variabili ausiliarie ζ = ka/2 e η = k 0 a/2. Le soluzioni si possono trovare graficamente dall’intersezione delle due curve: η2 + ζ 2 =

2m W a2 , ¯h2 4

η = ζ tan ζ

(2.57)

• Soluzioni dispari (α = π/2), per le quali k/ tan(ka/2) = −k 0 . Si procede come sopra, cercando le intersezioni delle curve η2 + ζ 2 =

2m W a2 , ¯h2 4

10

η = −ζ/ tan ζ

(2.58)

x tan x -x/tan x sqrt(1-x^2) sqrt(9-x^2) sqrt(36-x^2)

8

6

4

2

0 0

1.5708

3.1416

4.7124

6.2832

Si trova uno spettro (ovverosia l’insieme delle soluzioni) discreto, ovvero formato da valori isolati di E. C’e’ sempre almeno una soluzione (è una caratteristica del potenziale considerato che non vale per altre forme di potenziale). Lo spettro discreto è una caratteristica degli stati legati, ovvero confinati in una zona di spazio. Notiamo anche che: • Le soluzioni sono non degeneri; Questa è una caratteristica dello spettro discreto nei sistemi unidimensionali (teorema di non degenerazione). • Le soluzioni sono o pari: ψ(x) = ψ(−x), o dispari: ψ(x) = −ψ(−x), rispetto all’operazione di inversione, x → −x. Questa è una conseguenza della simmetria del potenziale, V (x) = V (−x). • La soluzione di pi` u bassa energia (lo stato fondamentale) è pari, quella di energia subito sopra è dispari, e cos`ı via. Questa è una proprietà generale dei potenziali unidimensionali simmetrici per inversione.

21

• Lo stato fondamentale non ha nodi (ovverosia non passa mai per lo zero: per nessun x, ψ(x) = 0); le soluzioni di energia crescente hanno un numero crescente di nodi (1,2,3,...). Anche questa è una proprietà generale dei potenziali unidimensionali. Ovviamente non dobbiamo dimenticare che esiste anche uno spettro continuo di soluzioni non legate e due volte degeneri per E > 0. E‘ utile il confronto con il caso della buca infinita, le cui soluzioni si ottengono banalmente imponendo che la funzione d’onda sia nulla cove il potenziale diventa infinito. Si tratta di condizioni diverse da quelle imposte per discontinuità finite del potenziale; ci se ne può convincere con una procedura di limite. Conviene traslare l’origine rispetto al caso precedente e consideraro un potenziale V (x) = 0 fra x = 0 e x = a, V (x) = ∞ al di fuori di tale intervallo. Imponiamo la condizione ψ(0) = ψ(a) = 0 sulle soluzioni per la particella libera. Si ottiene ψn (x) = sin(kn x),

kn =

nπ a

En =

¯h2 k 2 n 2 h2 = , 2m 8ma2

n = 1, ..., ∞ (2.59)

Si ritrovano tutte le caratteristiche del caso della barriera finita, salvo il numero di soluzioni (infinito per la barriera infinita) e lo spettro continuo (qui assente). Si notano inoltre due aspetti molto importanti: • le funzioni d’onda sono ortogonali fra di loro: Z

ψn∗ (x)ψm (x) = 0

se

n 6= m,

(2.60)

• le funzioni d’onda formano un insieme completo, ovvero qualunque funzione d’onda pu` o essere espressa come somma, in generale infinita, delle soluzioni dell’equazione di Schrödinger. Tali aspetti sono presenti anche nel caso della buca finita, in quanto derivano da propriet` a generali dell’equazioni di Schrödinger, ma sono in questo caso particolarmente visibili.

22

Chapter 3

L’equazione di Schr¨ odinger unidimensionale: soluzione analitica e numerica In questo capitolo verr` a descritta una metodologia per risolvere sia analiticamente che numericamente l’equazione di Schrödinger indipendente dal tempo in una dimensione per un oscillatore armonico; l’estensione dei metodi numerici ad altri tipi di potenziali non comporta particolari difficoltà.

3.1

L’oscillatore armonico

L’oscillatore armonico è uno dei problemi fondamentali della dinamica classica, e anche della meccanica quantistica. Rappresenta il pi` u semplice sistema modello in cui sono presenti delle forze attrattive, quindi è un importante riferimento per tutti i fenomeni vibrazionali. Ad esempio, le vibrazioni di un sistema di particelle fra loro interagenti può essere descritto, con una opportuna trasformazione di coordinate, in termini di modi normali di vibrazione, ciascuno dei quali è in pratica un oscillatore armonico indipendente dagli altri. Lo stesso vale in meccanica quantistica, dove per un sistema con stati legati le frequenze vibrazionali non sono altro (a meno della costante di Planck) le energie associate ai livelli energetici permessi. Attraverso lo studio quantistico dell’oscillatore armonico si possono quindi capire diverse cose relative alla quantizzazione, e alle funzioni d’onda degli stati legati. In questo capitolo esporremo i risultati principali della teoria dell’oscillatore armonico, e cercheremo di mostrare come impostare un codice di calcolo che permetta di risolvere numericamente la relativa equazione di Schrödinger. Il programma risultante potr` a poi facilmente essere modificato per inserire un potenziale di interazione diverso da quello quadratico, e permettere cos`ı lo studio di problemi che invece possono essere molto difficili da attaccare dal punto di vista analitico.

23

3.1.1

Unit` a

L’equazione di Schr¨ odinger di un oscillatore armonico unidimensionale è [utilizzando una notazione simile alla (2.4)] 1 d2 ψ 2m = − 2 E − Kx2 ψ(x) 2 dx 2 ¯h

(3.1)

dove m è la massa e K la costante di forza (la forza a cui è soggetta la massa è cioè F = −Kx, proporzionale allo spostamento e diretta verso l’origine). Classicamente a un tale oscillatore corrisponde una frequenza (angolare) s

ω=

K m

(3.2)

` conveniente passare ad unit` E a adimensionali (in cui lavorano i due programmi presentati in seguito): poniamo

ξ=

mK ¯h2

1/4

x

,

ε=

E ¯hω

(3.3)

[usando la definizione (3.2) per ω] ottenendo l’equazione equivalente d2 ψ ξ2 = −2 ε − dξ 2 2

!

ψ(ξ)

(3.4)

che è espressa in unit` a adimensionali.

3.1.2

Soluzione e livelli energetici

Come si pu` o facilmente verificare, per grandi ξ (tali da poter trascurare ε) le soluzioni della (3.4) devono avere l’andamento asintotico ψ(ξ) ∼ ξ n e±ξ

2 /2

(3.5)

dove n ha un qualsiasi valore finito. Il segno + nell’esponente deve però essere scartato a priori perchè darebbe luogo a soluzioni divergenti e quindi non normalizzabili (inoltre, l’intuizione stessa ci dice che la particella non dovrebbe tendere ad allontanarsi da ξ = 0, punto verso cui è diretta la forza). Sembra quindi conveniente provare a scorporare l’andamento asintotico desiderato ponendo 2 (3.6) ψ(ξ) = H(ξ)e−ξ /2 dove H(ξ) è una funzione che a grandi ξ si deve comportare in modo che 2 l’andamento sia determinato dal secondo fattore e−ξ /2 . H(ξ) non deve, in 2 particolare, crescere come eξ , altrimenti saremmo in presenza di una delle soluzioni che non desideriamo. Con la posizione (3.6) la (3.4) diventa, per la nuova funzione incognita H(ξ), H 00 (ξ) − 2ξH 0 (ξ) + (2ε − 1)H(ξ) = 0

24

(3.7)

Vediamo subito che ε0 = 1/2, H0 (ξ) = 1 è la soluzione pi` u semplice. Come tra poco si vedr` a, questa è la soluzione che rappresenta lo stato fondamentale, cioè quello ad energia pi` u bassa. Per ottenere le soluzioni generali sviluppiamo H(ξ) in una serie (in principio infinita): H(ξ) =

∞ X

An ξ n ,

(3.8)

n=0

deriviamo la serie per ottenere le derivate e riscriviamo la (3.7) combinando i termini con la stessa potenza di ξ: ∞ X

[(n + 2)(n + 1)An+2 + (2ε − 2n − 1)An ] ξ n = 0

(3.9)

n=0

Affinchè ci` o sia soddisfatto per qualsiasi valore di ξ è necessario che tutti i coefficienti siano nulli: (n + 2)(n + 1)An+2 + (2ε − 2n − 1)An = 0

(3.10)

Cos`ı, una volta dati A0 e A1 , la (3.10) permette di determinare per ricursione l’intera soluzione in forma di serie di potenze. Supponiamo che la serie sia veramente una serie infinita. A grandi n i termini si comportano quindi come An+2 2 → An n

(3.11)

2n Ma, ricordando che exp(ξ 2 ) = n ξ /n!, i cui coefficienti soddisfano pure alla (3.11), vediamo che questa relazione tra i coefficienti fa crescere H(ξ) come exp(ξ 2 ), ossia ci fornisce delle soluzioni divergenti indesiderate. L’unica maniera per evitare che questo accada è fare in modo che, nella (3.10), tutti i coefficienti da un certo punto in poi siano nulli, in modo che la serie si riduca in realtà ad un polinomio di grado finito. Questo avviene se e solo se

P

ε=n+

1 2

(3.12)

dove n è un intero positivo o nullo. Corrispondentemente, le energie possibili per l’oscillatore armonico sono quantizzate: 1 En = n + ¯hω n = 0, 1, 2, . . . (3.13) 2 I corrispondenti polinomi Hn (ξ) sono detti polinomi di Hermite. Hn (ξ) è di grado n in ξ, ha n nodi, ed è pari [Hn (−ξ) = Hn (ξ)] o dispari [Hn (−ξ) = 2 −Hn (ξ)] a seconda che n sia pari o dispari. Poichè e−ξ /2 non ha nodi ed è pari, anche l’intera autofunzione corrispondente all’autovalore dell’energia En ψn (ξ) = Hn (ξ)e−ξ

2 /2

(3.14)

ha n nodi e la parit` a di n. Pi` u sotto si mostra come una parità definita è una conseguenza della simmetria del problema rispetto all’inversione dell’asse x.

25

Figure 3.1: Funzioni d’onda e densità di probabilità dell’oscillatore armonico quantistico. I polinomi di Hermite di ordine pi` u basso sono H0 (ξ) = 1 ,

H1 (ξ) = 2ξ ,

H2 (ξ) = 4ξ 2 − 2 ,

H3 (ξ) = 8ξ 3 − 12ξ (3.15)

Un grafico delle corrispondenti funzioni d’onda e densità di probabilità è riportato in fig. 3.1.

3.1.3

Energia di punto zero

Una nota conseguenza della soluzione (3.13) è che il livello energetico pi` u basso—lo stato fondamentale—ha una energia finita ¯hω/2, chiamata energia di punto zero e tipica dei sistemi quantistici. La sua esistenza è legata al principio di indeterminazione di Heisenberg. Assumiamo infatti—in un’ottica semiclassica—che l’energia totale sia dell’ordine di (∆p)2 /2m+K(∆x)2 /2, dove ∆p e ∆x sono misure della dispersione tipica della quantità di moto e della posizione della particella. Il principio di indeterminazione ci dice che ∆x∆p ≥ ¯h/2, da cui possiamo estrarre √ ∆x ' ¯h/2∆p e minimizzare l’energia rispetto a ∆p. Si ottiene (∆p)2 ' ¯h Km/2, da cui E ' ¯hω/2. Dunque vediamo che l’energia minima non può essere nulla. Se lo fosse, avremmo determinato esattamente sia la posizione che la quantità di moto, in contraddizione col principio di indeterminazione. Le conseguenze dell’energia di punto zero possono essere importanti: ad esempio, He4 (a pressione atmosferica) resta allo stato liquido fino a temperature arbitrariamente piccole a causa dell’energia di punto zero.

3.1.4

Simmetria e parit` a

Tutte le autofunzioni dell’oscillatore armonico con n pari sono funzioni pari, e ` facile dimostrare che in casi come quelle con n dispari sono funzioni dispari. E questo in cui il potenziale è simmetrico, ossia V (−x) = V (x), una soluzione dell’equazione di Schr¨ odinger è necessariamente pari o dispari per motivi di simmetria.

26

Si immagini infatti di invertire l’asse x: x → −x. Nessuna osservabile fisica può cambiare per effetto di questa trasformazione, perchè il potenziale non varia. Poichè la densit` a di probabilità è un’osservabile, dovrà quindi essere |ψn (−x)|2 = |ψn (x)|

(3.16)

Ciò è possibile solo se le due funzioni differiscono per un fattore di fase complesso: ψn (−x) = eiα ψn (x) (3.17) con α reale. Effettuando due volte questa operazione di inversione dell’asse si ritorna però alla situazione di partenza. Quindi, applicando due volte in sequenza l’equazione qui sopra, si scopre che deve essere e2iα = 1

(3.18)

ossia α = mπ con m intero. La ψn è dunque pari se m è pari, e dispari se m è dispari. Si pu` o quindi a priori dire che, a causa della simmetria del potenziale, i polinomi di Hermite di grado pari devono avere tutti i coefficienti dispari nulli, e viceversa.

3.1.5

Confronto con la densit` a di probabilit` a classica

Le densit` a di probabilit` a delle funzioni d’onda ψn (x) dell’oscillatore armonico hanno, in generale, n + 1 picchi, la cui altezza aumenta mentre ci si avvicina ai corrispondenti punti di inversione classici. Queste densit` a di probabilit` a possono essere confrontate con quella dell’oscillatore armonico classico, in cui la massa si muove secondo x(t) = x0 sin(ωt). La probabilit` a ρ(x)dx di trovare la massa fra x e x + dx è proporzionale al tempo impiegato per attraversare quella regione, ossia inversamente proporzionale alla velocit` a, espressa in funzione di x: ρ(x)dx ∝

dx v(x)

(3.19)

q

Poichè v(t) = x0 ω cos(ωt) = ω x20 − x20 sin2 (ωt), sarà ρ(x) ∝ q

1 x20

(3.20)

− x2

Questa densit` a di probabilit` a ha un minimo a x = 0, e diverge ai punti di ` ovviamente nulla oltre il punto di inversione. inversione. E La densit` a di probabilit` a quantistica nello stato fondamentale è completamente diversa: presenta un massimo a x = 0, e decresce aumentando x. Al punto di inversione classico il suo valore è circa il 60% del valore massimo. La particella ha una elevata probabilità di trovarsi nella regione classicamente proibita. Nel limite di grandi numeri quantici, la densità quantistica tende tuttavia ad assomigliare a quella classica, ma esibisce il comportamento oscillatorio nella regione permessa tipico dei sistemi quantistici.

27

3.2 3.2.1

Meccanica quantistica e codici numerici: alcune considerazioni Energie cinetiche negative

Uno dei fatti nuovi importanti della meccanica quantistica rispetto a quella classica è la presenza di “energie cinetiche negative”, ossia la funzione d’onda può non essere nulla (e quindi la probabilità di trovare una particella essere finita) nelle “regioni proibite” dal punto di vista classico, V (x) > E. Basandoci sulla (2.4) e immaginando di essere in un caso semplice in cui V si può considerare costante, questo significa d2 ψ = k 2 ψ(x) (3.21) dx2 dove k 2 è un numero positivo, e questo a sua volta implica un comportamento esponenziale: sia ψ(x) ' exp(kx) che ψ(x) ' exp(−kx) soddisfano alla (3.21). Come sappiamo dagli studi delle buche di potenziale, generalmente solo una di queste due possibilit` a ha significato fisico: quella che dà luogo a una funzione d’onda che decresce esponenzialmente man mano che ci si addentra nella regione classicamente proibita. Queste regioni si traducono però spesso in serie difficoltà per i codici numerici, che per la loro natura generale contemplano entrambi i tipi di soluzioni. Tutti sappiamo che crescite esponenziali portano inevitabilmente a catastrofi, e cos`ı anche un algoritmo di integrazione tende a far esplodere la funzione d’onda in modo catastrofico nel momento in cui è presente, anche se in piccolissima quantit` a, una componente crescente nella soluzione; ed è inevitabile ` comune quindi che una funzione d’onda ottenuta numeriche questo accada. E camente, e perfettamente valida nella regione classicamente permessa, diverga improvvisamente se ci si addentra oltre un certo limite all’interno della regione classicamente proibita. In questo senso, trattare numericamente sistemi quantistici richiede pi` u attenzione che trattare sistemi classici, che sono intrinsecamente pi` u stabili.

3.2.2

Effetti della quantizzazione

Un importante punto da tener presente ai fini della risoluzione numerica di problemi quantistici, strettamente connesso al punto precedente, è la presenza di quantizzazione dei livelli energetici possibili per gli stati legati, espressa dall’espressione (3.13) nel caso dell’oscillatore armonico, ma fatto generale della meccanica quantistica. I livelli energetici possibili En non sono in generale noti a priori. Pertanto, in un’equazione di Schr¨ odinger (2.4) l’incognita non è solo ψ(x), ma anche E. Per ogni livello energetico, o autovalore En ci sarà una corrispondente funzione d’onda, o autofunzione ψn (x). Cosa succede se si cerca di risolvere un’equazione di Schrödinger utilizzando un’energia E non corrispondente ad un autovalore? La risposta che ci viene dallo studio delle autofunzioni dell’oscillatore armonico è che la quantizzazione delle energie nasce proprio imponendo le opportune condizioni al contorno, al

28

fine di impedire divergenze non fisiche della funzione d’onda nelle regioni proibite. Pertanto, se E non corrisponde ad un autovalore possiamo sicuramente aspettarci di osservare divergenze di ψ(x). I codici numerici che ricercano le energie permesse dovranno pertanto essere in grado di ”riconoscere” i problemi causati da energie sbagliate e saper aggiustare il tiro, modificando l’energia fino a portarla in coincidenza di un autovalore. Il programma presentato alla fine di questo capitolo implementa una strategia di questo genere.

3.3

Il metodo di Numerov

Vogliamo ora considerare il problema della risoluzione dell’equazione di Schrödinger indipendente dal tempo e unidimensionale in maniera numerica. Questo ci permetter` a di apprendere la metodologia generale da applicare per casi specifici (ad esempio quello dell’oscillatore armonico), e di comprendere la potenza e le limitazioni del metodo numerico. Il metodo di Numerov è utile per integrare equazioni differenziali del secondo ordine della forma generale d2 y = −g(x)y(x) + s(x) dx2 dove g(x) e s(x) sono funzioni date, e condizioni iniziali della forma y(x0 ) = y0

,

y 0 (x0 ) = y00

(3.22)

(3.23)

L’equazione di Schr¨ odinger (2.4) ha questa forma, con g(x) ≡ 2m [E − V (x)] h2 ¯ e s(x) = 0. Si vedr` a in seguito che anche le parti radiali di equazioni di Schrödinger in tre dimensioni a simmetria sferica appartengono a questa classe. Un’altra importante equazione che ricade in questa categoria è l’equazione di Poisson dell’elettromagnetismo, d2 φ = −4πρ(x) (3.24) dx2 dove ρ(x) è una densit` a di carica. In questo caso g(x) = 0 e s(x) = −4πρ(x). La metodologia è simile a quella dell’algoritmo di Størmer-Verlet per seguire l’evoluzione temporale di un punto materiale, sostituendo la coordinata spaziale a quella temporale. Dividiamo dunque l’intervallo spaziale di interesse in N intervallini di ampiezza ∆x, siano xi i nodi della griglia cos`ı ottenuta e yi = y(xi ) i valori della funzione incognita y(x) in corrispondenza di tali punti. Analogamente indichiamo con gi e si i valori delle funzioni (date) g(x) e s(x) negli stessi punti. Al fine di ottenere una equazione alle differenze finite sviluppiamo in serie di Taylor attorno ad un punto xn , spingendoci fino al quinto ordine: yn−1 = yn − yn0 ∆x + 12 yn00 (∆x)2 − 61 yn000 (∆x)3 + +O[(∆x)6 ] yn+1 = yn + yn0 ∆x + 12 yn00 (∆x)2 + 61 yn000 (∆x)3 + +O[(∆x)6 ]

1 0000 4 24 yn (∆x)

−

1 00000 5 120 yn (∆x)

1 0000 4 24 yn (∆x)

+

1 00000 5 120 yn (∆x)

(3.25)

29

Sommiamo le due equazioni: yn+1 + yn−1 = 2yn + yn00 (∆x)2 +

1 0000 y (∆x)4 + O[(∆x)6 ] 12 n

(3.26)

La (3.22) ci dice che yn00 = −gn yn + sn

(3.27)

Inoltre, indicando temporaneamente con zn questa quantità, sarà anche vero zn+1 + zn−1 = 2zn + zn00 (∆x)2 + O[(∆x)4 ]

(3.28)

zn+1 + zn−1 − 2zn + O[(∆x)2 ] (∆x)2

(3.29)

e quindi yn0000 ≡ zn00 =

Inserendo questi risultati nella (3.26) yn+1 = 2yn − yn−1 + (−gn yn + sn )(∆x)2 1 (−gn+1 yn+1 + sn+1 − gn−1 yn−1 + sn−1 + 2gn yn − 2sn )(∆x)2 + 12 +O[(∆x)6 ] (3.30) da cui la formula di Numerov h

yn+1 1 + gn+1 (∆x) 12

2

i

h

= 2yn 1 − 5gn (∆x) 12

2

i

h

− yn−1 1 + gn−1 (∆x) 12 2

6 +(sn+1 + 10sn + sn−1 ) (∆x) 12 + O[(∆x) ]

2

i

(3.31)

che permette di ottenere yn+1 a partire da yn e yn−1 e quindi ricorsivamente— dalle condizioni iniziali date—tutta la funzione sull’intervallo di interesse. Dalle condizioni iniziali (3.23) è ovviamente possibile integrare muovendosi sia nella direzione degli x positivi che in quella degli x negativi, e in presenza di simmetria rispetto ad un punto di inversione basterà integrare in una direzione sola. Nel caso di nostro interesse—l’equazione di Schrödinger—tutti i termini sn sono assenti, e in pratica è conveniente porre fn ≡ 1 + gn

(∆x)2 12

(3.32)

dove

2m [E − V (xn )] ¯h2 Con questa posizione la formula di Numerov si riduce a gn =

yn+1 =

(12 − 10fn )yn − fn−1 yn−1 fn+1

30

(3.33)

(3.34)

3.3.1

Programma: harmonic0

Il programma harmonic0.f901 (oppure harmonic0.c2 ) risolve l’equazione di Schrödinger per l’oscillatore armonico quantistico, integrando mediante l’algoritmo di Numerov descritto sopra, e ricercando gli autovalori mediante il metodo di ”shooting”. Si tratta di una procedura del tutto analoga a quella della ricerca dello zero di una funzione mediante il metodo di bisezione (descritto in C.1.1). Il programma va alla ricerca della soluzione ψn (x) con un numero di nodi assegnato n, e considera inizialmente l’energia E corrispondente al punto medio dell’intervallo [Emin , Emax ] (che siamo certi contenere l’autovalore desiderato En ). La funzione d’onda viene integrata partendo da x = 0 e muovendosi verso gli x positivi, e allo stesso tempo viene contato il numero di nodi, ossia di cambiamenti di segno della funzione. Se tale numero risulta essere superiore a n, significa che E è troppo alta; se invece il numero di nodi è minore o uguale a n, significa che E è troppo bassa. Viene allora scelto il semiintervallo opportuno—rispettivamente quello inferiore [Emin , E] o quello superiore [E, Emax ]—e la procedura iterata sul semiintervallo. Si considera di essere arrivati a convergenza quando l’ampiezza dell’intervallo di energia è scesa al di sotto di una soglia prefissata. Per x negativi la funzione viene costruita per simmetria, essendo ψn (−x) = (−1)n ψn (x). Questo è ovviamente possibile in quanto V (−x) = V (x), altrimenti sarebbe stato necessario integrare su tutto l’intervallo. Il programma chiede, nell’ordine, il valore massimo xmax a cui estendere l’integrazione (un valore tipico è compreso fra 5 e 10), il numero N di punti sulla griglia (da cui ∆x = xmax /N ), il nome del file di output e il numero di nodi richiesto. Ad un’ultima domanda sull’energia da provare si risponderà in generale con 0 per avviare il meccanismo di ricerca dell’autovalore descritto sopra; è tuttavia possibile anche inserire una specifica energia, in generale non corrispondente ad un autovalore, per forzare il programma ad effettuare una singola integrazione su quell’energia ed esaminare la funzione d’onda risultante. Questo permette di effettuare dei test per capire meglio il funzionamento del metodo di ricerca dell’autovalore. Il file di output contiene la soluzione finale ed è organizzato su cinque colonne contenenti rispettivamente x, ψ(x), |ψ(x)|2 , ρcl (x) e V (x). ρcl (x) è la densità di probabilit` a dell’oscillatore armonico classico data dalla (3.20) (normalizzata a 1). Oltre a questo file, il programma emette sullo standard output ad ogni iterazione il numero dell’iterazione, il numero di nodi trovati (sul solo semiasse degli x positivi) e la stima corrente dell’autovalore dell’energia. Come si vedr` a, è impossibile evitare che la soluzione diverga al di sopra di un certo x. Questo è il risultato dell’inevitabile presenza di una componente ∼ exp(+x2 /2) che, anche se piccola, porta ad una divergenza quando x sufficientemente grande. Tale divergenza comporta anche una difficoltà nel normalizzare la ψ(x). A causa di questa difficoltà, questo programma deve chiaramente essere migliorato. 1 2

http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/F90/harmonic0.f90 http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/C/harmonic0.c

31

3.3.2

Programma: harmonic1

Il programma harmonic1.f903 (oppure harmonic1.c4 ) costituisce una versione migliorata di harmonic0, che evita il problema della divergenza a grandi x. Il trucco utilizzato per evitare la divergenza è quello di effettuare due integrazioni: una in avanti partendo da x = 0, e una all’indietro partendo da xmax . L’autovalore viene fissato dalla condizione che il raccordo fra le due parti sia continuo con derivata prima continua. Il punto di raccordo scelto coincide col punto di inversione classica xcl , tale che V (xcl ) = E. Tale punto dunque varia man mano che diversi E vengono provati, ed è ovviamente determinato entro la precisione che può fornire la griglia. In pratica viene identificato l’indice icl del primo punto della griglia xc = icl∆x tale che V (xc ) > E. Possiamo quindi solo dire che il punto di inversione classico è situato fra xc − ∆x e xc . L’integrazione verso destra viene portata avanti fino a icl, ottenendo cos`ı una funzione ψL (x) definita in [0, xc ], contando ancora il numero di cambiamenti di segno n come in harmonic0. Se n non è quello corretto l’energia è abbastanza distante dal valore richiesto, e si procede ad un aggiustamento per bisezione come descritto in precedenza (senza che sia necessario integrare al di là di xc , regione in cui si sa a priori che non vi possono essere nodi). Se invece il numero di nodi è quello giusto, il programma procede con l’integrazione da destra verso sinistra5 , fermandosi allo stesso indice icl corrispondente a xc e ottenendo cos`ı una funzione ψR (x) definita in [xc , xmax ]. A quel punto vi sono due valori della funzione d’onda in xc : ψL (xc ) e ψR (xc ). La prima operazione effettuata è quella di riscalare ψR (x) di un fattore ψL (xc )/ψR (xc ), in modo che via sia un raccordo continuo fra le due funzioni in xc . Fatto questo, l’intera ψ(x) viene rinormalizzata in modo che sia R |ψ(x)|2 dx = 1. Ora inizia la parte nuova e cruciale: il calcolo della discontinuità della 0 (x ) − ψ 0 (x ). Questa differenza dovrebbe derivata prima alla giunzione, ψR c L c pure essere nulla per una buona soluzione, ma questo non sarà vero se non (all’interno di una precisione specificata) quando E = En . Il segno della differenza ci permette di capire se E è troppo alta o troppo bassa, e quindi di applicare nuovamente il metodo di bisezione per migliorarne la stima. Per calcolare la discontinuit` a, utilizziamo (indicando per brevità con i l’indice corrispondente al punto di griglia icl): L = y L − y 0L ∆x + 1 y 00L (∆x)2 + O[(∆x)3 ] yi−1 i i 2 i R = y R + y 0R ∆x + 1 y 00R (∆x)2 + O[(∆x)3 ] yi+1 i i 2 i

(3.35)

Notiamo che yiL = yiR = yi , e anche yi00L = yi00R = yi00 = −gi yi , come garantisce 3

http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/F90/harmonic1.f90 http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/C/harmonic1.c 5 Si noti lo statement y(mesh) = dx: ha il solo scopo di forzare il segno della soluzione ad essere positivo, in quanto la parte ottenuta a sinistra del punto di raccordo è pure positiva. Il valore dx non ha particolare importanza, dal momento che in seguito la soluzione viene riscalata per raccordarsi con l’altra in modo continuo alla giunzione. 4

32

il metodo di Numerov. Sommando fra loro le due equazioni abbiamo allora L R yi−1 + yi+1 = 2yi + (yi0R − yi0L )∆x − gi yi (∆x)2 + O[(∆x)3 ]

(3.36)

da cui

L + y R − [2 − g (∆x)2 ]y yi−1 i i i+1 + O[(∆x)2 ] ∆x o anche, utilizzando la notazione (3.32),

yi0R − yi0L =

yi0R − yi0L =

L + y R − (14 − 12f )y yi−1 i i i+1 + O[(∆x)2 ] ∆x

(3.37)

(3.38)

Il programma calcola in questo modo la discontinuità nella derivata prima. Se il suo segno è positivo, l’energia è troppo alta e viene quindi scelto il semiintervallo inferiore; se negativo, l’energia è troppo bassa e viene scelto il semiintervallo superiore. Al solito, si ritiene di essere giunti a convergenza quando l’ampiezza dell’intervallo di energie è diventata minore di una certa piccola tolleranza. Oltre alla funzione d’onda finale scritta sul file di output in modo analogo ad harmonic0 (ma ora è normalizzata!), il programma emette una linea sullo standard output per ogni iterazione, contenente quattro numeri: il numero dell’iterazione, il numero di nodi trovati (sul solo semiasse degli x positivi), la discontinuit` a nella derivata yi0R − yi0L (solo quando è stata calcolata, ossia se il numero di nodi era quello giusto—negli altri casi viene riportato zero) e la stima corrente dell’autovalore dell’energia.

3.3.3

Laboratorio

Ecco alcuni spunti per esperimenti numerici da effettuare in laboratorio, utilizzando entrambi i programmi: • Costruire e osservare6 le autofunzioni per diversi valori di n. • Esaminare le funzioni d’onda ottenute specificando una ben precisa energia non corrispondente ad un autovalore. Vedere la differenza di comportamento tra harmonic0 e harmonic1 in questi casi. • Osservare cosa succede quando l’energia è prossima ma non coincidente con un autovalore. Nuovamente, confrontare harmonic0 con harmonic1. • Esaminare gli effetti del ∆x (a quanti nodi si può arrivare?) e delle tolleranze. Alcuni spunti per modifiche del codice: • Modificare il potenziale, mantenendo la simmetria per inversione. 6

Un suggerimento: è utile graficare assieme alle autofunzioni o alle densit` a di probabilit` a anche la densit` a di probabilit` a classica nella quarta colonna, che indica la posizione dei punti di inversione classici. Con gnuplot questo si ottiene aggiungendo in coda al comando plot: ,"filename" using 1:4 with lines

33

• Modificare il potenziale, rompendo la simmetria per inversione. Si può per esempio considerare il potenziale seguente: h

i

V (x) = D e−2ax − 2e−ax + 1 .

(3.39)

E’ il potenziale di Morse, tuttora usato per modellizzare l’energia potenziale di una molecola biatomica. Che tipo di modifiche occorre apportare all’algoritmo per coprire questo caso? Esaminare lo spettro degli autovalori dell’energia e la forma delle funzioni d’onda. Confrontare con i risultati per l’oscillatore armonico con costante K scelta in modo tale da dare la stessa curvatura al minimo del potenziale di Morse: V (r) ' Da2 x2 ≡

k 2 x , 2

k = 2Da2 .

(3.40)

Verificare gli autovalori rispetto alla soluzione analitica esatta:

p

1 1 1 ¯hω − ¯hω n + 2 4D 2

E(n) = n +

2

(3.41)

p

dove ω = k/m = a 2D/m, o, in termini dei livelli Eh (n) dell’oscillatore armonico:

E(n) = Eh (n) 1 −

Eh (n) , 4D

Eh (n) = n +

1 ¯hω. 2

(3.42)

Potete per esempio considerare il caso D = 12.5, a = 0.2, in unità nelle quali 2m/¯h2 = 1.

34

Chapter 4

Propriet` a dell’equazione di Schr¨ odinger Questo capitolo descrive proprietà fondamentali dell’equazione di Schrödingere degli operatori associati a grandezze fisiche. Molti risultati sono ottenuti in modo pi` u formale nell’appendice A. Qui si usa invece il formalismo della funzione d’onda, pi` u diffuso e probabilmente pi` u semplice da capire.

4.1

Ortonormalit` a delle funzioni d’onda

Si dice che un insieme di funzioni d’onda complesse ψn (x) sono fra loro ortonormali se Z ψn∗ (x)ψm (x)dx = δnm (4.1) Nel caso n = m questo esprime semplicemente la normalizzazione a 1 di una funzione d’onda, mentre nel caso m 6= n esprime una condizione di ortogonalit` a1 fra le due funzioni. Dato l’operatore hamiltoniano H=−

¯h2 ∂ 2 + V (x), 2m ∂x2

(4.2)

le cui autofunzioni (funzioni d’onda soluzione dell’equazione di Schrödinger) hanno un valore determinato (autovalore) En dell’energia: Hψn (x) = En ψn (x),

(4.3)

vogliamo far vedere che le ψn (x) sono fra loro ortonormali, ossia soddisfano alla (4.1). Per dimostrarlo, si considerino le due equazioni Hψm = Em ψm H ∗ ψn∗ = En∗ ψn∗

(4.4)

La condizione di ortogonalit` a fra due vettori ordinari, ~a · ~b = i ai bi = 0 è analoga alla (4.1) qualora si consideri l’integrale come una “somma su tutti gli x”, e che la coniugazione complessa del primo termine è necessaria per far s`ı che il prodotto scalare di un oggetto con se stesso sia una misura della sua “norma”. 1

P

35

Si moltiplichi a sinistra la prima per ψn∗ e la seconda per ψm e si integri: ψn∗ Hψm dx = EmR ψn∗ ψm dx ψm H ∗ ψn∗ dx = En∗ ψm ψn∗ dx

R R

R

(4.5)

Gli integrali nei membri a destra sono identici. Facciamo ora vedere che i membri a sinistra devono essere uguali fra loro. Ciò corrisponde a dimostrare l’hermiticit` a dell’operatore H, ovvero la seguente proprietà: Z

ψn∗ Hψm dx =

Z

∗ ψm Hψn dx

∗

.

(4.6)

L’identit` a è palese per quanto riguarda la parte di H relativa al potenziale, perchè V (x) è un semplice fattore moltiplicativo reale. Concentriamoci allora sul termine cinetico ed effettuiamo una integrazione per parti: Z

x2 d2 ψm ∗ dψm dx = ψ − n dx2 dx x1

Z

d2 ψn∗ dψn∗ x2 dx = ψ − ψm m dx2 dx x1

Z

ψn∗

dψn∗ dψm dx dx dx

(4.7)

dψm dψn∗ dx dx dx

(4.8)

e analogamente Z

Assumiamo che agli estremi di integrazione x1 e x2 la funzione d’onda e le sue derivate sia nulla (se cos`ı non fosse, si immagini di racchiudere l’intero sistema, per quanto grande sia, in una scatola limitata da barriere infinite, e di far corrispondere gli estremi di integrazione con queste barriere). Se i termini integrati sono nulli, abbiamo trovato che Z

d2 ψm dx = dx2

Z

ψn∗ Hψm dx =

Z

ψn∗

d2 ψn∗ dx dx2

(4.9)

ψm H ∗ ψn∗ dx

(4.10)

ψm

e quindi che Z

La (4.5) ci d` a allora (Em − En∗ )

Z

ψn∗ ψm dx = 0

(4.11)

Quando m = n, l’integrale è 1, e quindi deve essere En∗ = En

(4.12)

ossia gli autovalori dell’energia sono senz’altro reali (da notare che tutte le grandezze fisiche misurabili, in quanto reali, sono rappresentate da operatori hermitiani). Supponiamo che sia m 6= n. Se anche Em 6= En , allora deve essere Z

ψn∗ ψm dx = 0

(4.13)

ossia autofunzioni corrispondenti ad autovalori diversi sono sicuramente ortogonali fra loro. Pu` o per` o anche accadere che sia m 6= n ma Em = En : è il caso della degenerazione.

36

In questo caso la (4.13) potrebbe non essere soddisfatta, ma si possono sempre scegliere le autofunzioni in modo che lo sia. Supponiamo ad esempio che le prime p autofunzioni ψ1 , . . . , ψp appartengano allo stesso autovalore E e siano linearmente indipendenti (ossia nessuna di esse può essere espressa come combinazione lineare delle altre), ma non siano fra loro ortogonali. Ogni loro combinazione lineare è quindi anch’essa una soluzione appartenente allo stesso autovalore E. A partire dalle ψ1 , . . . , ψp possiamo allora costruire un nuovo insieme di autofunzioni del tutto equivalente ψˆ1 , . . . , ψˆp , in cui le autofunzioni sono ortogonali fra loro. Si pu` o ad esempio procedere in questo modo: ψˆ1 (x) = ψ1 (x) ψˆ2 (x) = a1 ψˆ1 (x) + a2 ψ2 (x) ψˆ3 (x) = b1 ψˆ1 (x) + b2 ψˆ2 (x) + b3 ψ3 (x) ...

(4.14)

dove a1 e a2 sono determinati da Z

0=

ψˆ1∗ ψˆ2 dx = a1

Z

|ψˆ1 |2 dx + a2

Z

ψˆ1∗ ψ2 (x)dx

(4.15)

e dalla condizione di normalizzazione, e cos`ı via (si tratta del metodo di ortonormalizzazione di Gram-Schmidt). Possiamo perci` o interpretare la (4.13) nel senso che se gli autovalori delle due autofunzioni sono diversi esse devono essere ortogonali, mentre se sono uguali esse possono essere scelte in modo da essere ortogonali.

4.2

Sviluppo di una soluzione generica

La linearit` a dell’equazione di Schrödinger ci assicura che, se ψn (x) sono le autofunzioni dell’equazione indipendente dal tempo, tutte le loro combinazioni lineari X Ψ(x, t) = cn ψn (x)e−iEn t/¯h (4.16) n

sono anche una soluzione. La (4.1) ci permette di dimostrare l’inverso, ossia che una generica soluzione Ψ(x, t) della (2.8) può sempre essere sviluppata secondo la (4.16) (propriet` a di completezza). Per dimostrare ci` o, consideriamo l’istante t = 0 e poniamo Z

cm =

∗ ψm (x)Ψ(x, 0)dx

(4.17)

e quindi consideriamo la funzione X

cn ψn (x)

(4.18)

n

Questa funzione è una soluzione (in quanto combinazione lineare di soluzioni), e all’istante t = 0 deve coincidere con Ψ(x, t). Infatti vale per ogni m l’identità cm =

X n

cn δmn =

X

Z

cn

∗ ψm (x)ψn (x)dx

Z

=

" ∗ ψm (x)

# X n

n

37

cn ψn (x) dx (4.19)

P

[si noti che abbiamo usato la (4.1)]. Questo è possibile solo se n cn ψn (x) coincide con Ψ(x, 0). Ma dato il valore di una funzione d’onda ad un certo istante, la sua evoluzione è completamente determinata dall’equazione di Schrödinger P temporale, e quindi l’evoluzione di Ψ(x, 0) dovrà essere la stessa di n cn ψn (x), cioè quella indicata dalla (4.16).

4.3

Valori medi

Uno dei postulati della meccanica quantistica è che il valore medio di una qualsiasi grandezza fisica funzione delle coordinate generalizzate qi e dei corrispondenti momenti generalizzati (quantità di moto) pi , F (q, p), si ottiene costruendo un operatore differenziale dove ad ogni pi nell’espressione classica si sostituisce −i¯h∂/∂qi . Il valore medio di questa grandezza su uno stato descritto dalla funzione d’onda Ψ(q, t) sar` a allora dato da hF i =

Z

Ψ∗ F Ψdq

(4.20)

dove F (un operatore hermitiano) agisce sulla funzione che sta alla sua destra, e l’integrazione è effettuata su tutte le coordinate del sistema. Se F è solo funzione delle coordinate, questa regola non è sorprendente: ci R dice che il valore medio è F (q)|Ψ|2 dq, ossia una normale media pesata sulla densit` a di probabilit` a. La situazione pi` u interessante riguarda le dipendenze dalle quantit` a di moto. Consideriamo ad esempio l’energia cinetica T per una particella in tre dimensioni: T (p) = p2 /2m. Applicando la regola, otteniamo hT i = −

¯h2 2m

dove ∇2 =

Z

Ψ∗ ∇2 Ψ dx dy dz

∂2 ∂2 ∂2 + 2+ 2 2 ∂x ∂y ∂z

(4.21)

(4.22)

Nel caso di una particella libera, l’energia cinetica è l’intero hamiltoniano e quindi l’equazione di Schr¨ odinger è −

¯h2 2 ∂Ψ ∇ Ψ(r, t) = i¯h 2m ∂t

(4.23)

con soluzione (a meno di un fattore di normalizzazione) Ψ(r, t) = ei(k·r−Et/¯h)

(4.24)

dove E è una costante, e k un vettore costante legato ad E da ¯h2 k 2 =E 2m

(4.25)

La (4.24) rappresenta un’onda piana che si propaga lungo la direzione del vettore d’onda k, con lunghezza d’onda λ = 2π/k e frequenza ω = E/¯h.

38

Applichiamo a questo caso il metodo per trovare la quantità di moto dell’onda, ad esempio lungo la direzione x: hpx i = −i¯h

Z

Ψ∗ (r, t)

∂ Ψ(r, t) dx dy dz = h ¯ kx ∂x

(4.26)

In realt` a una soluzione del genere ha un valore ben definito, ossia è una autofunzione, dell’operatore quantità di moto. Il punto fondamentale è che mentre l’informazione relativa alle coordinate si trova nelle ampiezze delle funzioni d’onda, l’informazione relativa alla quantità di moto si trova nelle fasi. Gli operatori differenziali in pratica accedono alla fase.

4.4

La formulazione matriciale

Abbiamo visto che ogni soluzione dell’equazione di Schrödinger può essere espressa nella forma (4.16). Note le autofunzioni dell’hamiltoniano, una funzione d’onda è pertanto caratterizzata dai coefficienti (complessi) c0 , c1 , c2 , . . .

(4.27)

Possiamo pensare a questi numeri come alle componenti di un vettore in uno spazio dove ogni asse rappresenta una autofunzione dell’energia. Poichè la Ψ deve essere normalizzata, occorre che sia X

|cn |2 = 1

(4.28)

n

Qual è l’energia media di Ψ? Utilizzando la (4.20) per l’operatore H, tenendo conto di (4.3) e (4.1), si trova hEi =

X

|cn |2 En

(4.29)

n

Questo va interpretato nel senso che una misura dell’energia della funzione darà sempre come risultato uno degli autovalori En , e a ciascuno di essi è associata una probabilit` a |cn |2 . Supponiamo ora invece di voler calcolare il valor medio associato ad un generico operatore F : hF i = = =

Z

Ψ∗ (x, t)F Ψ(x, t)dx

Z X

c∗n ψn∗ (x)eiEn t/¯h F

(4.30) X

cm ψm (x)e−iEm t/¯h dx

n m X ∗ i(En −Em )t/¯ h cn cm Fnm e nm

dove si è definito Fnm ≡

Z

ψn∗ (x)F ψm (x)dx

(4.31) (4.32)

(4.33)

Questa quantit` a è detta elemento di matrice dell’operatore F fra le autofunzioni ψn e ψm . Cos`ı, in questa formulazione (sviluppata da Heisenberg prima della

39

scoperta della meccanica ondulatoria) gli stati sono rappresentati da vettori, e gli operatori da matrici. In un’altra notazione—sviluppata da Dirac—questo elemento di matrice è indicato con < n|F |m >

(4.34)

dove |m > rappresenta un autostato (generalizzazione del concetto di autofunzione) dell’hamiltoniano. La (4.32) ci dice che un valor medio, in generale, dipende dal tempo. Vediamo anche che se Fnm fosse una matrice diagonale il valor medio sarebbe una costante, in quanto non vi sarebbe pi` u alcun fattore di fase oscillante con coefficiente non nullo. I fattori di fase oscillanti con frequenza (En − Em )/¯h danno dei termini che sono legati a transizioni del sistema da uno stato ad un altro. Termini di questo tipo si trovano ad esempio quando si studiano i processi di emissione o ` da notare che—grazie assorbimento di radiazione elettromagnetica (fotoni) E alla doppia sommatoria su m e su n—tutti i termini sono in realtà reali in quanto per ogni termine viene sommato anche il suo complesso coniugato, corrispondente allo scambio di indici. Questo consente di continuare a interpretare i coefficienti che moltiplicano gli Fnm come delle probabilità. Un operatore F pu` o essere applicato ad una funzione. In questa rappresentazione matriciale, questa operazione corrisponde ad applicare una matrice a P un vettore, ottenendo un altro vettore. Infatti, se ψ = n cn ψn (e ricordando la regola (4.17) per ottenere le “componenti del vettore”): Z

(F ψ)m = Z

=

∗ ψm F ψdx

(4.35)

∗ ψm F

cn ψn dx

(4.36)

∗ ψm F ψn dx

(4.37)

X n

=

X

Z

cn

n

=

X

Fmn cn

(4.38)

n

che è la consueta regola dell’algebra lineare. Notiamo anche che si pu` o sviluppare F ψm =

X

a` ψ`

(4.39)

`

dove i coefficienti sono ottenuti con la consueta regola: Z

a` =

ψ`∗ (F ψm )dx = F`m

(4.40)

e quindi F ψm =

X `

40

F`m ψ`

(4.41)

Possiamo infine applicare due operatori F e G in sequenza, e mostrare come nella rappresentazione matriciale questa operazione corrisponda ad effettuare un prodotto tra le due matrici corrispondenti secondo le consuete regole dell’algebra lineare. Infatti, usando la (4.41) due volte: Z

ψn∗ (x)F Gψm (x)dx

(F G)nm = =

XZ

ψn∗ (x)F G`m ψ` (x)dx

(4.42) (4.43)

`

=

X

=

X

Z

G`m

ψn∗ (x)F ψ` (x)dx

(4.44)

ψn∗ (x)Fk` ψk (x)dx

(4.45)

`

Z

G`m

k`

=

X

Fk` G`m δkn

(4.46)

Fn` G`m

(4.47)

k`

=

X `

che è appunto l’ordinaria regola per il prodotto di matrici.

4.5

Regole di commutazione

Il valore medio del prodotto di operatori dipende dall’ordine in cui gli operatori vengono applicati. Prendiamo ad esempio una coordinata x e la quantità di moto ad essa coniugata p, e calcoliamo il valor medio del prodotto px nella rappresentazione delle coordinate: ∂ hpxi = Ψ −i¯h (xΨ) dx ∂x Z ∂Ψ ∗ = −i¯h Ψ Ψ + x dx ∂x = −i¯h + hxpi Z

∗

(4.48) (4.49) (4.50)

ovvero hxp − pxi = h[x, p]i = i¯h

(4.51)

dove con la notazione [A, B] indichiamo l’operatore AB − BA che chiameremo commutatore tra A e B. Diremo che A e B commutano quando il loro commutatore è nullo, ossia quando è indifferente l’ordine con cui sono applicati su uno stato. Come abbiamo appena visto questo non è sempre vero. Il risultato (4.51) non dipende dallo stato, ed è quindi una identità a livello di operatore: [x, px ] = i¯h (4.52) (dove si è aggiunto l’indice x in px per sottolineare che si tratta della quantità di moto coniugata a x). Si pu` o far vedere che le variabili come x e px che non commutano sono quelle non misurabili simultaneamente. Invece, [x, y] = 0

41

(4.53)

(non vi sono vincoli alla determinazione simultanea di diverse coordinate di posizione), [px , py ] = 0 (4.54) (lo stesso per le quantit` a di moto), e [x, py ] = 0

(4.55)

(lo stesso per la coordinata in una direzione e la quantità di moto in un’altra). In generale, date due quantità osservabili A e B rappresentate in meccanica quantistica da operatori, |h[A, B]i|/2 rappresenta il limite inferiore al prodotto ∆A∆B, dove ∆A e ∆B sono gli scarti quadratici medi di misure effettuate simultaneamente su queste due variabili. Si tratta di una versione pi` u generale del principio di indeterminazione. Se A e B non commutano, è impossibile determinarle entrambe simultaneamente con precisione assoluta. D’altra parte, si pu` o vedere che se Φ è un’autofunzione comune di A e B: AΦ = aΦ

,

BΦ = bΦ

(4.56)

(dove a e b sono gli autovalori, ossia dei semplici numeri) allora ABΦ = AbΦ = bAΦ = baΦ

(4.57)

BAΦ = BaΦ = aBΦ = abΦ

(4.58)

e sono uguali, ossia [A, B]Φ = 0. Se questo è vero per un insieme completo di autofunzioni (ad esempio, per tutte le autofunzioni dell’energia ψn ), allora ne segue necessariamente [A, B] = 0. Si può dimostrare che è vero anche l’inverso: se A e B commutano, allora hanno un insieme completo di autofunzioni in comune.

4.6

Quantit` a conservate

Vogliamo ora dimostrare che dato un operatore F , e definito l’operatore dF/dt in modo tale che per ogni stato dipendente dal tempo Ψ si abbia

dF dt

=

d hF i dt

(4.59)

ossia il valor medio di dF/dt sullo stato sia pari alla derivata temporale del valor medio di F sullo stesso stato, allora vale la importante relazione dF = [F, H] (4.60) dt Questa relazione ci permette di identificare facilmente le quantit` a conservate, che cioè non variano nel tempo: sono quelle che commutano con l’hamiltoniano. Per dimostrarlo, consideriamo la (4.32) per il valor medio di F su uno stato Ψ sviluppato come somma di autofunzioni dell’energia, e deriviamola rispetto al tempo, ottenendo cos`ı i¯h

dF dt

=

X d i hF i = c∗n cm (En − Em )Fnm ei(En −Em )t/¯h dt ¯ h nm

42

(4.61)

Questo ci consente di identificare gli elementi di matrice dell’operatore dF/dt:

dF dt

= nm

i (En − Em )Fnm ¯h

(4.62)

(questo fa s`ı che la (4.32) valga anche per questo operatore, come deve essere!). Costruiamo ora invece l’elemento di matrice dell’operatore [F, H]: [F, H]nm = (F H − HF )nm Z

= Z

=

ψn∗ F Hψm dx −

(4.63) Z

ψn∗ F Em ψm dx −

ψn∗ HF ψm dx

(4.64) !

Z

ψn∗ H

X

F`m ψ` dx

(4.65)

`

Z

= Em

ψn∗ F ψm dx −

X

Z

F`m

ψn∗ E` ψ` dx

(4.66)

`

= Em Fnm −

X

F`m E` δn`

(4.67)

`

= (Em − En )Fnm

(4.68)

Confrontando la (4.62) con la (4.68) vediamo che

i¯h

dF dt

= [F, H]nm

(4.69)

nm

Ma se questo vale per tutti gli elementi di matrice, l’uguaglianza deve avvenire a livello di operatore, ossia la (4.60) deve essere vera. Siamo ora in posizione di comprendere cosa dobbiamo fare per classificare in un modo utile gli stati di un sistema quantistico. • La prima cosa da fare è cercare gli autovalori En dell’hamiltoniana H. Questo ci fornisce dei numeri quantici n utili ai fini della classificazione. Sappiamo che le autofunzioni di H corrispondono a sono stazionari [vedere (2.13)], quindi questi numeri quantici non variano nel tempo: sono buoni numeri quantici. • Possono per` o esserci delle degenerazioni: ad un certo En possono corrispondere diversi stati. Questi stati differiranno per altri numeri quantici, che vorrei saper determinare. • Devo allora cercare un altro operatore A che commuti con H: [A, H] = 0. Questo garantisce che i suoi autovalori siano costanti nel tempo; e anche che siano determinabili con esattezza e simultaneamente agli autovalori di H. • Un solo operatore addizionale potrebbe non bastare a classificare gli stati. Cercher` o allora un altro operatore B, che deve soddisfare anch’esso a [B, H] = 0. Ma non basta! Occorre anche che sia [A, B] = 0. Se cos`ı non fosse, non potrei determinare simultaneamente gli autovalori di A e di B, e quindi un tale schema sarebbe inutile ai fini della classificazione.

43

• Ripeto il procedimento finchè sono riuscito a classificare tutti gli stati, e non esistono altri operatori che commutino con H e tutti gli altri. • Ho allora costruito un insieme di osservabili che determina univocamente lo stato del sistema. Nel caso dell’atomo di idrogeno, come si vedrà, si utilizzano quattro operatori: H, L2 , Lz e Sz , discussi in seguito.

44

Chapter 5

Atomi con un elettrone 5.1

Equazione di Schr¨ odinger in un campo centrale

Consideriamo un sistema quantistico costituito da due particelle di masse m1 e m2 interagenti tra loro, e in assenza di campi esterni. Supponiamo per il momento che il potenziale di interazione V (r) sia arbitrario, anche se sappiamo che nel caso dell’atomo di idrogeno l’interazione è coulombiana. Vogliamo trovare prima i risultati generali del problema che non dipendono dalla natura specifica del potenziale. Il potenziale V non pu` o comunque dipendere che dalla sola distanza |r2 −r1 | tra le due particelle, e l’hamiltoniano sarà H=

p2 p21 + 2 + V (|r2 − r1 |) 2m1 2m2

(5.1)

Come in meccanica classica, si può effettuare un cambiamento di variabili e passare alle due nuove variabili m1 r1 + m2 r2 m1 + m2 r = r2 − r1

R =

(5.2) (5.3)

` corrispondenti alla posizione del centro di massa e alla posizione relativa. E conveniente anche definire M

= m1 + m2 m1 m2 m = m1 + m2

(5.4) (5.5)

dove m è detta massa ridotta. Si pu` o facilmente vedere che, definendo anche i nuovi operatori corrispondenti P = −i¯h∇R e p = −i¯h∇r , l’hamiltoniano diventa H=

P2 p2 + + V (r) 2M 2m

(5.6)

da cui si vede immediatamente che le variabili si separano. Il moto del centro di massa è quello di una particella libera di massa M ; la soluzione è un’onda piana.

45

La parte interessante è ovviamente quella relativa. L’equazione di Schrödinger corrispondente è la stessa che avrebbe una massa m immersa in un campo di forze centrali V (r), con simmetria sferica rispetto all’origine. Nel caso degli atomi con un elettrone, l’interazione è fra il protone (o un nucleo pi` u pesante) e l’elettrone, e quindi il rapporto fra le masse è pari ad almeno 1836. La massa ridotta sar` a quindi appena pi` u piccola di quella dell’elettrone. L’equazione di Schr¨ odinger che studieremo in questo capitolo è allora: #

"

¯h2 2 ∇ + V (r) ψ(r) = Eψ(r) Hψ(r) ≡ − 2m

5.2

(5.7)

Il momento angolare

La soluzione classica del problema di una particella in un campo centrale (ossia soggetta ad un potenziale V (r) dipendente solo dalla distanza rispetto a un punto fisso) passa attraverso l’introduzione di una quantità, il momento angolare (o momento della quantit` a di moto), definita come L=r×p

(5.8)

dove r è il vettore posizione e p il vettore quantità di moto. In meccanica classica si trova che L è una quantit` a conservata, con importanti conseguenze tra cui la planarit` a dell’orbita. Ci aspettiamo che anche il corrispondente operatore quantistico giochi un ruolo importante, ed infatti cos`ı è. Possiamo immediatamente dire qualcosa sulle sue proprietà di commutazione, facendo uso delle (4.52)—(4.55) e utilizzando la proprietà generale (immediatamente dimostrata) [AB, C] = A[B, C] + [A, C]B

(5.9)

Si trova [Lx , x] = 0

,

[Lx , y] = i¯hz

,

[Lx , z] = −i¯hy

(5.10)

[Lx , py ] = i¯hpz

,

[Lx , pz ] = −i¯hpy

(5.11)

e [Lx , px ] = 0

,

e propriet` a analoghe ottenute ciclando gli indici per Ly e Lz . Si può far vedere che analoghe propriet` a valgono per i commutatori fra componenti di L: [Lx , Lx ] = 0

,

[Lx , Ly ] = i¯hLz

,

[Lx , Lz ] = −i¯hLy

(5.12)

e in realt` a è vero per qualsiasi grandezza vettoriale A, funzione arbitraria di coordinate e quantit` a di moto: [Lx , Ax ] = 0

,

[Lx , Ay ] = i¯hAz

,

[Lx , Az ] = −i¯hAy

(5.13)

Inoltre, dati due vettori A e B (sempre corrispondenti ad operatori quantistici), si può costruire l’operatore “prodotto scalare” A · B = Ax Bx + Ay By + Az Bz

46

(5.14)

e risulta [Lx , A · B] = [Ly , A · B] = [Lz , A · B] = 0

(5.15)

come si dimostra subito usando le (5.13). In particolare, facendo coincidere A e B con L stesso, abbiamo anche [Lx , L2 ] = [Ly , L2 ] = [Lz , L2 ] = 0

(5.16)

Come si vedr` a nella sezione 5.4, e come intuibile dal risultato classico, per una particella in un campo centrale L2 commuta con H, ed è quindi una quantità conservata che d` a origine a un buon numero quantico. Anche ogni singola componente di L commuta con H. Però, le (5.12) mostrano che due diverse componenti di L non commutano fra loro, e non sono pertanto misurabili simultaneamente.

5.3

Autofunzioni del momento angolare

Esprimiamo il momento angolare nella rappresentazione delle coordinate: L = −i¯hr × ∇

(5.17)

Consideriamo un sistema di riferimento polare (r, θ, φ), dove l’asse polare coincide con l’asse cartesiano z, θ è l’angolo polare e φ quello azimutale. Siano ur , uθ e uφ i versori (che costituiscono una terna ortonormale destrorsa) associati a spostamenti in cui varia solo r, θ o φ rispettivamente. Si ha ∇ = ur

∂ 1 ∂ 1 ∂ + uθ + uφ ∂r r ∂θ r sin θ ∂φ

(5.18)

Applicando la (5.17), 1 ∂ ∂ L = −i¯h uφ − uθ ∂θ sin θ ∂φ

(5.19)

Esprimendo i versori della terna polare in funzione di quelli della terna cartesiana ur = sin θ cos φ ux + sin θ sin φ uy + cos θ uz

(5.20)

uθ = cos θ cos φ ux + cos θ sin φ uy − sin θ uz

(5.21)

uφ = − sin φ ux + cos φ uy

(5.22)

possiamo calcolare le componenti cartesiane di L nello spazio polare. In particolare risulta ∂ Lz = −i¯h (5.23) ∂φ e " # 1 ∂ ∂ 1 ∂2 2 2 L = −¯h sin θ + (5.24) sin θ ∂θ ∂θ sin2 θ ∂φ2

47

Cerchiamo ora le autofunzioni dell’operatore L2 , che torneranno utili in seguito risolvendo l’equazione di Schrödinger per una particella in un campo centrale: (5.25) L2 Y (θ, φ) = ¯h2 `(` + 1)Y (θ, φ) dove abbiamo espresso in questo modo (per futura convenienza) l’autovalore. Notiamo che, moltiplicando i due membri per − sin2 θ/¯h2 , l’equazione agli autovalori diventa ∂Y (θ, φ) ∂ sin θ sin θ ∂θ ∂θ

+

∂ 2 Y (θ, φ) = −`(` + 1) sin2 θ Y (θ, φ) ∂φ2

(5.26)

Supponiamo che la soluzione sia separabile in una funzione di solo θ e una di solo φ: Y (θ, φ) = Θ(θ)Φ(φ) (5.27) e dividiamo il risultato per ΘΦ: 1 ∂ ∂Θ sin θ sin θ Θ ∂θ ∂θ

+ `(` + 1) sin2 θ = −

1 ∂2Φ Φ ∂φ2

(5.28)

Il primo membro è funzione solo di θ, e il secondo solo di φ. Entrambi devono allora essere uguali a una costante positiva1 , che indichiamo con m2 . Abbiamo allora ottenuto due equazioni: 1 d dΘ sin θ sin θ dθ dθ

#

"

m2 Θ=0 + `(` + 1) − sin2 θ

d2 Φ + m2 Φ = 0 dφ2

(5.29)

(5.30)

La seconda ci dice che deve essere Φ(φ) = Ce±imφ

(5.31)

Poichè φ è un angolo azimutale, è necessario che m sia intero affinchè la funzione sia ad un solo valore. La (5.29), usando cos θ come variabile, è nota in fisica matematica come equazione di Legendre. Si pu` o risolvere in modo analogo a quanto fatto per l’oscillatore armonico: esprimendo cioè la soluzione in forma di una serie di potenze di cos θ, e richiedendo che non diverga per alcun valore di cos θ. Risulta che una divergenza a cos θ = 1 può essere evitata solo se si assume che la serie sia in realt` a un polinomio di grado finito, ossia che tutti i coefficienti da un certo grado in poi siano nulli. Si può vedere che questo implica ` intero, e ` ≥ |m|. Le funzioni risultanti sono indicate con P`m (cos θ) e si chiamano funzioni associate ai polinomi di Legendre. I polinomi di Legendre P` (cos θ) sono le soluzioni dell’equazione di Legendre per m = 0, e le funzioni associate sono ad essi connesse da

P`m (w) = 1 − w2 1

m/2 dm

dwm

P` (w).

(5.32)

Si pu` o verificare che un valore negativo porterebbe a soluzioni esponenziali non accettabili

48

Le autofunzioni dell’operatore L2 hanno dunque la forma Y`m (θ, φ) = C`m P`m (cos θ)eimφ

(5.33)

dove C`m è una costante di normalizzazione, e sono dette armoniche sferiche. Poiche Lz , dato dalla (5.23), opera solo su φ, queste sono anche autofunzioni di questo operatore: (5.34) Lz Y`m (θ, φ) = ¯hmY`m (θ, φ) In sostanza, h ¯ 2 `(` + 1) rappresenta il modulo quadrato del momento angolare, e ¯hm la sua proiezione lungo l’asse z. ` dev’essere un intero positivo o nullo, e m un intero compreso fra −` e `. Per un dato ` ci sono dunque 2` + 1 valori permessi per m.

5.4

Separazione in parte radiale e angolare

Torniamo al nostro problema di una particella in un campo centrale. Introduciamo, analogamente a quanto fatto nella sezione 5.2, un sistema di riferimento polare (r, θ, φ), dove l’operatore gradiente è dato dalla (5.18), e l’operatore laplaciano (come si pu` o far vedere con un po’ di pazienza) da ∂ 1 ∂ r2 ∇ = 2 r ∂r ∂r 2

1 ∂ ∂ + 2 sin θ r sin θ ∂θ ∂θ

+

1 ∂2 r2 sin2 θ ∂φ2

(5.35)

Confrontando con la (5.24), si vede che questo si può scrivere 1 ∂ ∂ ∇ = 2 r2 r ∂r ∂r

2

−

L2 r2 ¯h2

(5.36)

dove L2 dato dalla (5.24) contiene esclusivamente termini dipendenti dagli angoli. Possiamo allora scrivere l’hamiltoniano come ¯h2 1 ∂ ∂ r2 H=− 2m r2 ∂r ∂r

+

L2 + V (r) 2mr2

(5.37)

Un termine L2 /2mr2 appare anche nell’analogo problema classico: altri non è che il “potenziale centrifugo”, ossia un potenziale fittizio che genera una “forza” che tende ad allontanare la massa dall’origine, e che discende dal fatto che il sistema è in rotazione (se L2 > 0), e che noi stiamo osservando la sola variabile radiale. Classicamente si può dunque tener conto dell’effetto della rotazione considerando un potenziale efficace Vˆ (r) = V (r) + L2 /2mr2 , dove il secondo termine tende a spingere la massa verso gli r crescenti. Vediamo ora la situazione nel caso quantistico. Un’ispezione della forma (5.37) ci mostra subito che [L2 , H] = 0 (5.38) che ci garantisce che L2 è conservato, ossia i suoi autovalori non dipendono dal tempo e i due operatori hanno un insieme di autofunzioni in comune. Già

49

sappiamo quindi che gli autovalori di L2 potranno essere usati per classificare gli stati. L’espressione (5.23) per Lz e la regola (5.16) ci dicono anche che [Lz , H] = 0

(5.39)

e quindi anche gli autovalori di Lz saranno conservati e potranno essere usati per classificare gli stati. Procediamo ora alla separazione della variabile radiale da quelle angolari, la cui possibilit` a è fortemente suggerita sia da questi risultati che dall’osservazione della (5.37). Poniamo ψ(r, θ, φ) = R(r)Y (θ, φ) (5.40) quindi riscriviamo l’equazione di Schrödinger (5.7), dividendola per RY : 1 ∂ ∂R ¯h2 r2 − 2 2m r R(r) ∂r ∂r

+

1 L2 Y + V (r) = E 2mr2 Y

(5.41)

o ancora moltiplicando per −2mr2 /¯h2 e riarrangiando, 1 ∂ ∂R r2 R(r) ∂r ∂r

−

1 2 2mr2 2 [V (r) − E] = 2 L Y ¯h ¯h Y

(5.42)

Il membro sinistro dipende solo da r, quello destro solo da θ e φ, e quindi entrambi devono essere uguali ad una costante. Abbiamo già [vedi (5.25)] indicato questa costante con `(` + 1), e trovato che ` deve essere un intero affinchè la soluzione non diverga. Le soluzioni per la parte angolare sono le armoniche sferiche Y`m (θ, φ) date dalla (5.33). Dovrà quindi essere 1 ∂ ∂R r2 R(r) ∂r ∂r

−

2mr2 [V (r) − E] = `(` + 1) ¯h2

(5.43)

ovvero l’equazione di Schr¨ odinger per la parte radiale è ¯h2 1 ∂ ∂Rn` − r2 2m r2 ∂r ∂r

#

"

¯h2 `(` + 1) + V (r) + Rn` (r) = En` Rn` (r) 2mr2

(5.44)

Ci aspettiamo che in generale le energie dipendano da ` perchè il potenziale efficace dipende da `; inoltre per un dato ` ci aspettiamo per gli stati legati (se ve ne sono!) una quantizzazione dei livelli energetici, e abbiamo indicato con n il corrispondente indice. La funzione d’onda totale sarà allora ψn`m (r, θ, φ) = Rn` (r)Y`m (θ, φ)

(5.45)

L’energia non dipende da m. Come già osservato, m rappresenta la proiezione del momento angolare su un asse scelto arbitrariamente. A causa della simmetria sferica del problema, l’energia non può dipendere dall’orientamento del vettore L, ma solo dal suo modulo. All’energia En` sarà dunque associata una degenerazione 2` + 1 (o maggiore, se esistono altri osservabili commutanti che non abbiamo considerato!).

50

5.4.1

Funzioni d’onda angolari

Le funzioni d’onda angolari per un problema a simmetria sferica non dipendono dunque dalla natura del potenziale, e sono date dalle armoniche sferiche Y`m (θ, φ) (5.33). Il loro aspetto per diversi valori di ` e m può essere esaminato ad esempio nella “galleria” dell’università di Oviedo2 , oppure esplorato attivamente usando l’applet Java al Davidson College3 . Si noti che m rappresenta la proiezione del momento angolare sull’asse z. Pertanto, le funzioni con m = 0 tenderanno a essere disposte lungo tale asse; quelle con m = ` tenderanno a localizzarsi prevalentemente sul piano xy. Le armoniche sferiche di ordine pi` u basso sono le seguenti: Y00 (θ, φ) = Y11 (θ, φ) = Y10 (θ, φ) = Y22 (θ, φ) = Y21 (θ, φ) = Y20 (θ, φ) =

q q q

1/4π

(5.46)

3/8π sin θ eiφ

(5.47)

3/4π cos θ

(5.48)

q

15/32π sin2 θ e2iφ

(5.49)

q

15/8π sin θ cos θ eiφ

q

5/16π 3 cos2 θ − 1 .

(5.50) (5.51)

Si sono assunte funzioni normalizzate secondo la normalizzazione tradizionale: Z

∗ Ylm (θ, φ)Ylm (θ, φ)dΩ = δll0 δmm0

(5.52)

dove Ω è l’angolo solido. L’ortogonalità delle armoniche sferiche è una naturale conseguenza del loro carattere di autofunzioni del momento angolare (nonché dell’equazione che soddisfano). Considerare −m al posto di m significa cambiare il segno all’esponente del termine exp(imφ) ossia, in pratica, a prendere la funzione complessa coniugata. E’ però da notare che la fase delle armoniche sferiche è arbitraria e che esistono diverse convenzioni Per identificare il valore di ` viene spesso usata la notazione spettroscopica: si indicano con s, p, d, f , g, . . . rispettivamente gli stati con ` = 0, 1, 2, 3, 4, . . .

5.5

Il potenziale coulombiano

Il caso pi` u importante e famoso è quello in cui V (r) è il potenziale coulombiano: V (r) = −

Ze2 , 4π0 r

(5.53)

dove e = 1.6021 × 10−19 C è la carica dell’elettrone, Z è il numero atomico (numero di protoni nel nucleo), 0 = 8.854187817 × 10−12 in unità MKSA. In 2 3

http://www.unioviedo.es/qcg/harmonics/harmonics.html http://webphysics.davidson.edu/Applets/Hydrogenic/

51

fisica si usa ancora molto il sistema CGS, nel quale il potenziale coulombiano ha la forma:: V (r) = −Zqe2 /r. (5.54) u semplice Nel seguito si user` a qe2 = e2 /(4π0 ) in modo da ricondursi alla pi` forma CGS. ` spesso comodo lavorare in unit` E a atomiche (a.u.): le unità di lunghezza sono espresse in raggi di Bohr (o semplicemente bohr), a0 : a0 =

¯h2 = 0.529177˚ A = 0.529177 × 10−10 m, me qe2

(5.55)

mentre le energie sono espresse in Rydberg (Ry): 1Ry =

me qe4 = 13.6058eV. ¯h2

(5.56)

dove me è la massa dell’elettrone, non la massa ridotta. E’ immediato verificare che in tali unit` a, ¯h = 1, me = 1/2, qe2 = 2. Se invece del Ry si prende l’Hartree (Ha): 1 Ha = 2 Ry =

me qe4 = 27.212 eV ¯h2

(5.57)

come unit` a di energia, si ottiene un altro set di unità atomiche, nelle quali ¯h = 1, me = 1, qe = 1. Attenzione alla confusione! Mai parlare di ”unità atomiche” senza specificare chiaramente quali. Nel seguito si useranno occasionalmente le prime (unit` a atomiche ”Rydberg”).

5.6

La funzione d’onda radiale per atomi idrogenoidi

` conveniente porre E χ(r) = rR(r)

(5.58)

e scrivere l’equazione radiale per χ(r) anzichè R(r). Si vede facilmente che la (5.44) diventa "

#

¯h2 d2 χ Zqe2 ¯h2 `(` + 1) + E + − χ(r) = 0 2me dr2 r 2me r2

(5.59)

Notiamo come questa equazione sia del tutto analoga all’equazione di Schrödinger in una dimensione (2.4), per una particella soggetta ad un potenziale efficace Zq 2 ¯h2 `(` + 1) Vˆ (r) = − e + (5.60) r 2me r2 Come gi` a sottolineato, il secondo termine è il potenziale centrifugo. Gli stessi metodi utilizzati per trovare la soluzione della (2.4) (e in particolare il metodo numerico di Numerov) possono quindi essere utilizzati per trovare le autofunzioni radiali dell’energia.

52

Notiamo innanzitutto che per piccoli r il potenziale centrifugo è il termine dominante del potenziale. L’andamento delle soluzioni per r → 0 sarà allora determinato da `(` + 1) d2 χ ' χ(r) (5.61) dr2 r2 che d` a χ(r) ∼ r`+1 , oppure χ(r) ∼ r−` . La seconda possibilità va scartata, perchè χ(r) non pu` o divergere. Per grandi r invece, notiamo che avremo stati legati se E < 0 (in quanto esister` a un punto di inversione classico al di là del quale l’energia cinetica diventa negativa, e quindi la funzione d’onda decade esponenzialmente, e quindi solo alcune energie potranno dare luogo a soluzioni valide), e liberi se E > 0. Il caso E > 0 corrisponde a un problema di scattering elettrone-nucleo con uno spettro continuo di energie, e non ce ne occupiamo. L’andamento delle soluzioni per r → ∞ sar` a allora determinato da 2me d2 χ ' − 2 Eχ(r) (5.62) 2 dr ¯h √ che d` a χ(r) ∼ exp(±kr), dove k = −2me E/¯h. Il segno + va però scartato perchè comporta una divergenza indesiderata. Sembra allora sensato assumere per la soluzione una forma χ(r) = r`+1 e−kr

∞ X

An r n

(5.63)

n=0

che garantisce un comportamento corretto in entrambi i casi limite, purchè la serie non diverga esponenzialmente. L’equazione per l’atomo idrogenoide può essere risolta seguendo lo stesso procedimento utilizzato per l’oscillatore armonico nella sezione 3.1. Ossia, si inserisce lo sviluppo (5.63) nella (5.59), si trova una formula di ricorrenza per i coefficienti An , si fa vedere che la serie in generale diverge come exp(2kr) a meno che non si interrompa dando origine a un polinomio, e si fa infine vedere che questo accade solo in corrispondenza a particolari valori di E. In particolare questo accade per Z 2 me qe4 Z2 Ry En = − 2 = − (5.64) n 2¯h2 n2 dove n ≥ ` + 1 è un intero detto numero quantico principale. Per un dato ` si avranno quindi soluzioni per n = ` + 1, ` + 2, . . .; oppure, pensando fissato n, i valori possibili per ` sono ` = 0, 1, . . . , n − 1. La soluzione per la funzione d’onda radiale si scrive s

χn` (r) =

(n − ` − 1)!Z `+1 −x/2 2`+1 x e Ln+1 (x) n2 [(n + `)!]3 a30

(5.65)

dove si è posto s

2Z r 2me En x≡ =2 − r n a0 ¯h2

53

(5.66)

dove gli L2`+1 n+1 (x) sono i polinomi di Laguerre, di grado n − ` − 1. Il coefficiente è stato scelto in modo da ortonormalizzare l’insieme di funzioni: Z ∞ 0

χn` (r)χn0 ` (r)dr = δnn0

(5.67)

Abbiamo gi` a dimostrato che l’ortogonalità è garantita per le autofunzioni di un hamiltoniano a cui corrispondono autovalori diversi dell’energia [vedi (4.13)]. Sottolineiamo alcuni risultati rilevanti:

5.6.1

Densit` a radiale

La probabilit` a di trovare la particella a una distanza compresa tra r e r + dr dal centro è ottenuta integrando sulle variabili angolari: Z

dr

|ψn`m (r, θ, φ)|2 rdθ r sin θ dφ = |Rn` |2 r2 dr = |χn` |2 dr

(5.68)

avendo sfruttato la propriet` a di normalizzazione delle armoniche sferiche Z

|Y`m (θ, φ)|2 dθ sin θ dφ = 1

(5.69)

(dove l’integrazione è estesa a tutti i possibili angoli). Ne segue anche che la condizione di normalizzazione in termini di χ è Z ∞ 0

|χn` (r)|2 dr = 1

(5.70)

La funzione |χ(r)|2 pu` o essere dunque direttamente interpretata come una densit` a radiale.

5.6.2

Stato fondamentale

Lo stato fondamentale ha n = 1, e quindi ` = 0. Si tratta dunque del caso in cui il momento angolare è nullo, e la corrispondente armonica sferica è costante: l’autofunzione è quindi a simmetria sferica. L’energia dello stato per l’atomo di idrogeno (Z = 1) è pari a −1, ossia l’energia di legame dell’elettrone è pari ad un Rydberg (a parte la piccola correzione legata alla massa ridotta). La funzione d’onda dello stato fondamentale è, con la normalizzazione esatta, un semplice esponenziale: Z 3/2 ψ100 (r, θ, φ) = √ e−Zr/a0 π

5.6.3

(5.71)

Comportamento vicino al nucleo

Il termine dominante vicino al nucleo è quello corrispondente al primo termine della serie, ossia χn` (r) ∼ r`+1 . Quindi pi` u ` è grande, pi` u rapidamente la funzione tende a zero avvicinandosi al nucleo. Questo corrisponde al fatto che la funzione è “spinta via” dal potenziale centrifugo. Quindi le funzioni radiali con grande ` non penetrano vicino al nucleo.

54

5.6.4

Comportamento lontano dal nucleo

A grandi valori di r il comportamento è dominato dall’ultimo termine della serie, ossia va come χ(r) ∼ rn exp(−Zr/na0 ). Questo significa che (trascurando gli altri termini) |χn` (r)|2 ha un massimo attorno a r = n2 a0 /Z. Questo fornisce una stima grossolana della “dimensione” dell’autofunzione. La dimensione globale è dunque determinata soprattutto da n.

5.6.5

Numero di nodi

Poichè nella (5.65) compare un polinomio di grado n − ` − 1, questo è anche il numero di nodi della funzione. In particolare, le autofunzioni con ` = 0 hanno n − 1 nodi; e tutte quelle con ` = n − 1 non hanno nodi. L’aspetto delle funzioni radiali pu` o essere esaminato ad esempio sul sito di Wolfram Research4 o esplorato attraverso l’eccellente applet Java al Davidson College5 .

5.7

Degenerazione accidentale e simmetria dinamica

Nonostante il potenziale efficace che appare nella (5.59) dipenda da `, e la parte angolare delle autofunzioni pure dipenda assai fortemente da `, l’espressione (5.64) dipende solo da n. Abbiamo dunque una degenerazione delle energie sugli n possibili valori per `, che si aggiunge a quella di ordine 2` + 1 legata ai possibili valori del numero quantico m [implicata dalla (5.44) in cui m non appare]. La degenerazione complessiva6 associata a n è n−1 X

(2` + 1) = n2

(5.72)

`=0

La degenerazione delle energie per diversi valori di ` è una situazione molto particolare che si verifica soltanto quando il potenziale di interazione è coulombiano. Si tratta di cioè di una degenerazione accidentale, che scompare appena il potenziale non è pi` u puramente coulombiano. Una degenerazione indica generalmente la presenza di una simmetria, e quindi di una quantit` a conservata. Ad esempio la degenerazione in m è legata alla simmetria sferica e alla conservazione del momento angolare. Si può far vedere che il corrispondente classico della degenerazione accidentale negli atomi idrogenoidi è la conservazione del vettore di Runge-Lenz M=

p×L α − r m r

(5.73)

verificata per una hamiltoniana classica H=

p2 α − 2m r

4

http://library.wolfram.com/webMathematica/MSP/Explore/Physics/Hydrogen http://webphysics.davidson.edu/Applets/Hydrogenic/ 6 Come si vedr` a in seguito, in realt` a c’è ancora un fattore 2 dovuto allo spin. 5

55

(5.74)

` questo il caso del moto relativo di due corpi attratti dalla forza gravitazionale. E Come ben noto, le orbite sono ellittiche, e sono orbite chiuse: l’orientazione dell’ellisse non cambia nel tempo. Il vettore di Runge-Lenz è diretto lungo l’asse maggiore dell’ellisse. Il corrispondente vettore quantistico ha una espressione lievemente diversa ma sostanzialmente simile: M=

1 Zqe2 (p × L − L × p) − r 2m r

(5.75)

e si pu` o far vedere che M è ortogonale a L, e [M, H] = 0: ossia è una quantità conservata.

5.8

Programma: hydrogen

Il programma hydrogen radial.f907 oppure hydrogen radial.c8 risolve l’equazione ` sostanzialmente basato su harmonic1, con radiale per un atomo idrogenoide. E piccole differenze dovute all’equazione leggermente diversa, e alla risoluzione su griglia a passo logaritmico.

5.8.1

Griglia logaritmica

La risoluzione numerica diretta della (5.59) presenta qualche difficoltà a causa della singolarit` a del potenziale a r = 0. Le difficoltà si possono aggirare lavorando su una griglia a passo variabile in r anzichè costante, che diventa sempre pi` u fitta man mano che ci si avvicina all’origine. Una descrizione pi` u approfondita dello schema qui presentato si può trovare nel cap.6 di questo libro: C. Froese Fischer, The Hartree-Fock method for atoms, Wiley, 1977. Chiamiamo x la nuova variabile di integrazione. Definiremo una griglia a passo costante in x, in modo da poter continuare ad adottare il metodo di Numerov senza modifiche. In generale per una mappatura definita da x = f (r)

(5.76)

∆x = f 0 (r)∆r

(5.77)

f (r) ≡ log Zr

(5.78)

avremo La nostra particolare scelta è

che fornisce quindi ∆r (5.79) r Il rapporto ∆r/r si mantiene pertanto costante sulla griglia logaritmica cos`ı definita. ∆x =

7 8

http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/F90/hydrogen radial.f90 http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/C/hydrogen radial.c

56

Tuttavia, trasformando la (5.59) nella variabile x appare anche un termine con la derivata prima, che impedisce l’applicazione del metodo di Numerov. Il problema si aggira trasformando anche la funzione incognita in questo modo: 1 y(x) = √ χ (r(x)) r

(5.80)

` facile vedere che trasformando la (5.59) in modo da esprimerla in funzione di x E e y (ma il potenziale pu` o essere lasciato funzione di r), i termini con la derivata prima si cancellano, e moltiplicando l’intera equazione per r3/2 si ottiene "

Zqe2 2me 2 d2 y r + E + dx2 r ¯h2

!

1 − `+ 2

2 #

y(x) = 0

(5.81)

Come si vede, questa equazione non presenta pi` u singolarità per r = 0. Dal punto di vista pratico, la subroutine do mesh definisce all’inizio e una √ volta per tutte i valori di r, r, r2 per ogni punto della griglia. Il potenziale è pure calcolato una volta per tutte da init pot. La griglia è calcolata a partire da un x minimo pari a -8, corrispondente a Zrmin ' 3.4 × 10−3 unità atomiche (raggi di Bohr).

5.8.2

Applicazione della teoria perturbativa

Particolare attenzione merita la sezione del programma: i = icl ycusp = (y(i-1)*f(i-1)+f(i+1)*y(i+1)+10.d0*f(i)*y(i)) / 12.d0 dfcusp = f(i)*(y(i)/ycusp - 1.d0) * Aggiornamento autovalore usando teoria delle perturbazioni de = dfcusp/ddx12 * ycusp*ycusp * dx il cui scopo è quello di stimare al primo ordine in teoria delle perturbazioni la differenza δe fra l’autovalore attuale e quello esatto. Vogliamo spiegare il meccanismo che sta alla base di questa stima (prima di affrontare questa sezione occorre essere familiari col metodo perturbativo descritto nel prossimo capitolo). Ricordiamo che icl è l’indice corrispondente al punto di inversione classico. L’integrazione viene effettuata in avanti sino a questo indice, ed all’indietro pure sino a questo indice. icl è quindi l’indice di raccordo tra le due funzioni. La funzione di destra viene riscalata in modo che non vi sia una discontinuit` a della funzione y al punto di raccordo, tuttavia la derivata prima dy/dx sar` a in generale discontinua (a meno che l’autovalore non sia quello cercato). Si pu` o quindi dire che y(icl) è il valore previsto dalla ricorrenza di Numerov calcolata usando come punto centrale sia icl-1 che icl+1. La ricorrenza di Numerov usando come punto centrale icl non è però mai stata applicata. Il valore calcolato per ycusp è appunto quanto predice la ricorrenza di Numerov usando icl come punto centrale; il problema è che è diverso da y(icl). L’ottica in cui ci porremo ora è quella di pensare che la funzione ottenuta sia la soluzione esatta di un problema diverso, e precisamente di un problema in

57

cui il potenziale nel solo intervallo ∆x centrato su y(icl) è diverso da quello dato. Il potenziale modificato “curva” la soluzione nel modo osservato. Una volta trovata la modifica del potenziale necessaria per ottenere questo effetto, la teoria delle perturbazioni ci fornisce una stima della differenza di autovalore rispetto all’originale. Riesaminiamo la formula di Numerov (3.34), e notiamo come la formula (pensiamola applicata al caso n=icl-1) ci fornisce in realtà il solo prodotto y(icl)f(icl). Solitamente da questo prodotto estraiamo y(icl) perchè f(icl) si pensa dato. Supponiamo ora invece che f al punto icl abbia un diverso valore fcusp non noto, ma tale che la funzione trovata soddisfi alla relazione di Numerov anche nel punto icl. Deve allora essere fcusp ycusp = f(icl) y(icl) in quanto questo prodotto è quanto fornito dal metodo di Numerov (integrando da icl-1 o icl+1), e ycusp è il valore che y deve avere affinchè la ricorrenza di Numerov sia soddisfatta anche in icl. Questo spiega la definizione di dfcusp (variazione di fcusp) calcolata dal programma. Il prossimo step è quello di calcolare la corrispondente variazione del potenziale. Per “potenziale” intendiamo qui il potenziale unidimensionale efficace che appare nell’equazione (5.81), W (x) = (2me /¯h2 )r2 V + (` + 1/2)2 , e da qui nelle definizioni per gli f(i). Differenziando si trova che la variazione del potenziale δW che d` a luogo a una variazione δf è data da δW = −

¯h2 12 δf 2me (∆x)2

(5.82)

La teoria delle perturbazioni al primo ordine ci dà allora la corrispondente variazione dell’autovalore: δe = hy|δW |yi =

Z

|y(x)|2 δW (x)dx = |y(xc )|2 δW ∆x

(5.83)

avendo sviluppato l’integrale in una sommatoria sulla griglia, dove l’unico contributo non nullo viene fornito dalla regione di ampiezza ∆x centrata sul punto di inversione classico xc (ossia, nel programma, l’indice icl). Abbiamo dunque trovato che ¯h2 12 δe = − |y(xc )|2 δf ∆x (5.84) 2me (∆x)2 è la differenza di autovalore per passare dal potenziale originale a quello perturbato. A parte il diverso segno—in quanto noi desideriamo risalire all’autovalore del problema non perturbato—questa è esattamente l’espressione usata dal programma per il calcolo della correzione all’autovalore de. Finchè si è lontani dalla convergenza, questa correzione è grossolana e può portare a errori considerevoli. Per questo motivo il programma è “protetto” dalla linea e = max(min(e+de,eup),elw) che impedisce in tutti i casi che venga stimato un nuovo valore di e all’esterno dell’intervallo [elw,eip]. Quando il programma procede verso la convergenza, la stima diventa via via sempre migliore e consente una convergenza assai rapida nella fase finale.

58

5.8.3

Laboratorio

• Esaminare le soluzioni al variare di n e `; verificare la presenza della degenerazione accidentale. • Provare a modificare leggermente il potenziale definito nella subroutine init pot rimuovendo il suo carattere coulombiano, e verificare che la degenerazione accidentale scompare. Alcune modifiche possibili sono V (r) = −Zqe2 /r1+δ dove δ è un numero piccolo, positivo o negativo; oppure l’aggiunta di uno smorzamento esponenziale: V (r) = −Zqe2 /r exp(−Qr) dove Q è un numero dell’ordine di 0.05 a.u..

59

Chapter 6

Metodi approssimati La risoluzione esatta, analitica, dell’equazione di Schrödinger è possibile solo in pochi casi, e anche la risoluzione numerica, specie per sistemi a molte particelle, può rivelarsi non banale o non fattibile in pratica. Esistono tuttavia dei metodi approssimati estremamente utili, con cui si può in molti casi ridurre il problema completo a un caso pi` u semplice.

6.1

Metodo perturbativo

Supponiamo che l’operatore hamiltoniano possa essere scritto in forma H = H0 + V

(6.1)

dove H0 è un “operatore hamiltoniano di riferimento” le cui autofunzioni si suppongono note: 0 0 0 H0 ψ m = Em ψm (6.2) e V una “perturbazione” che supporremo piccola, e indipendente dal tempo. 0 siano tutti Supporremo inoltre che non vi sia degenerazione, ossia che gli Em diversi (il caso in cui vi è degenerazione si può pure trattare senza troppa difficolt` a). Consideriamo le soluzioni del problema completo: Hψn = En ψn

(6.3)

0 costituiscono un insieme completo, ` Poichè le ψm e di certo possibile sviluppare le ψn secondo X 0 ψn = anm ψm (6.4) m

Con questa posizione la (6.3) diventa X

X

0 anm En ψm

(6.5)

0 0 anm (En − Em )ψm

(6.6)

0 anm (H0 + V )ψm =

m

m

da cui X m

0 anm V ψm =

X m

60

Moltiplichiamo a sinistra per ψk0∗ e integriamo. Utilizziamo per brevità la notazione di Dirac Z 0 0 hψk0 |V |ψm (6.7) i ≡ ψk0∗ ψm dv (integrato su tutte le coordinate). Si ottiene, sfruttando l’ortonormalità delle 0 , ψm X 0 (6.8) anm hψk0 |V |ψm i = ank (En − Ek0 ) m

Fino a questo momento non è stata effettuata alcuna approssimazione. Poniamo ora En = En0 + n

(6.9)

anm = δnm + αnm

(6.10)

ottenendo hψk0 |V |ψn0 i +

X

0 αnm hψk0 |V |ψm i = (δnk + αnk )(En0 − Ek0 + n ) (6.11)

m

= n δnk + αnk (En0 − Ek0 + n ) (6.12) e supponiamo che n e αnm , assieme agli integrali (6.7), siano delle quantità piccole. Approssimiamo al primo ordine (eliminando cioè tutti i termini di ordine successivo al primo): hψk0 |V |ψn0 i = n δnk + αnk (En0 − Ek0 )

(6.13)

Nel caso k = n abbiamo n = hψn0 |V |ψn0 i

(6.14)

Questo è un risultato molto importante, che ci permette di ottenere rapidamente una stima del cambiamento di autovalore conseguente a una perturbazione mediante integrali sugli stati imperturbati. Il caso k 6= n ci dà αnk =

hψk0 |V |ψn0 i En0 − Ek0

(6.15)

che pure ha una interpretazione fisica: uno stato ψn è ottenuto come “miscela” degli stati imperturbati, di cui ψn0 costituisce l’ingrediente fondamentale, e in cui gli altri ψk0 appaiono con coefficienti proporzionali all’”accoppiamento” tra ψn0 e ψk0 attraverso la perturbazione, e inversamente proporzionali alla differenza di energia. Notiamo che al primo ordine deve essere αnn = 0. Infatti 1 = hψn |ψn i =

X

|δnm + αnm |2 = |1 + αnn |2 +

m

X

|αnm |2

(6.16)

m6=n

Al primo ordine il membro a destra è 1 + 2αnn , da cui l’asserto. Questo equivale ad affermare che la correzione alla funzione d’onda al primo ordine, ψn1 , è ortogonali alla funzione d’onda imperturbata, ψn0 .

61

Riassumendo, i risultati della teoria delle perturbazioni non dipendenti dal tempo al primo ordine sono (per autovalori non degeneri): En = En0 + hψn0 |V |ψn0 i ≡ En0 + En1 X hψ 0 |V |ψ 0 i m n 0 ψm ≡ ψn0 + ψn1 ψn = ψn0 + 0 − E0 E m n m6=n

(6.17) (6.18)

` naturalmente possibile sviluppare la teoria delle perturbazioni ad ordini E pi` u elevati. Ad esempio al secondo ordine si trova En = En0 + hψn0 |V |ψn0 i +

X |hψ 0 |V |ψ 0 i|2 m n . m6=n

0 En0 − Em

(6.19)

L’osservatore attento noter` a come tale espressione sia equivalente alla seguente: En = En0 + hψn0 |V |ψn0 i + hψn0 |V |ψn1 i,

(6.20)

ovvero: per avere la correzione al prim’ordine dell’energia, bastano le funzioni d’onda imperturbate; per la correzione al second’ordine dell’energia, servono le correzioni al prim’ordine delle funzioni d’onda. In generale, si dimostra il teorema 2n+1: la correzione all’ordine 2n+1 dell’energia si calcola dalla correzione fino all’ordine n sulle funzioni d’onda.

6.1.1

Perturbazioni con autovalori degeneri

Nel caso in cui uno o pi` u autovalori siano degeneri, la 6.15 ci mostra chiaramente che tale approccio ha un problema. In effetti le formule sopra ricavate si applicano solo se vale la condizione 0 0 |hψn0 |V |ψm i| 0) e vettore d’onda k = ω/c, possiamo prendere la seguente forma per il potenziale vettore: A i(k·r−ωt) (6.28) ee + e∗ e−i(k·r−ωt) 2 dove e è il vettore di norma unitaria ortogonale a k, in generale complesso, che rappresenta la polarizzazione della radiazione e.m.. Introduciamo ora l’Eq. (6.28) nell’Eq.(6.24). Poniamo ωf i = (Ef − Ei )/¯h,

A(r, t) =

eA hEf |(pn ·e)eik·rn |Ei i 2mc¯h e osserviamo il risultato:

Gn =

Fn =

eA hEf |(pn ·e)e−ik·rn |Ei i (6.29) 2mc¯h

2 Z t X Z t i(ωf i −ω)t0 0 i(ωf i +ω)t0 0 Pf i (t) = Fn e dt + Gn e dt . 0 0

(6.30)

n

Avremo termini oscillanti, la cui media temporale va a zero, salvo che quando ¯hω = Ef − Ei (per il primo termine) o h ¯ ω = Ei − Ef (per il secondo termine). Nel primo caso abbiamo assorbimento di radiazione e Ef > Ei ; nel secondo caso abbiamo emissione stimolata, con Ef < Ei . Il calcolo pu` o essere portato a termine e la probabilit` a di transizione per unit` a di tempo determinata: 2

Wf i

X dPf i (t) 4π 2 q 2 ≡ = 2 e2 I(ωf i ) hEf | (pn · e)eik·rn |Ei i , dt m c¯hωf i n

Ei < Ef ,

(6.31) formula valida sia per l’assorbimento i → f che per l’emissione stimolata f → i con intensit` a della luce I. Il calcolo dell’emissione spontanea - che non appare nel quadro della teoria semiclassica - si può fare con considerazioni sull’equilibrio nel corpo nero e porta ad una formula in cui appare di nuovo lo stesso elemento di matrice al quadrato che nella (6.31), con coefficienti a moltiplicare differenti.

6.2.1

Transizioni di dipolo

Nella sezione precedente abbiamo ritrovato le condizioni note sull’energia della luce incidente affinché una transizione possa avvenire. Tuttavia tale condizione è solo necessaria. Resta da determinare il valore dell’elemento di matrice: X e Df i = hEf | (pn · e)eik·rn |Ei i (6.32) m n che appare nella (6.31), che se nullo (o molto piccolo) renderà la transizione proibita. Per frequenze tipiche di transizione elettroniche negli atomi, k ∼ 1/1000˚ A, mentre r ∼ 1˚ A. Di conseguenza, k · rn 0; • sulla proiezione del momento angolare: definiamo gli operatori D+ = Dx + iDy , D− = Dx − iDy . Gli elementi di matrice non nulli sono: hEf , L0 , M 0 |D+ |Ei , L, M i se M 0 = M + 1; hEf , L0 , M 0 |D− |Ei , L, M i se M 0 = M − 1; hEf , L0 , M 0 |Dz |Ei , L, M i se M 0 = M . Quest’ultima regola è particolarmente utile perché permette di distinguere le varie transizioni fra stati con M differente, usando luce polarizzata: infatti lo specifico elemento di matrice del dipolo dipende dalla polarizzazione e.

6.3

Metodo variazionale

Consideriamo un operatore hamiltoniano H e una funzione ψ, che può essere fatta variare liberamente con la condizione che resti normalizzata a 1. Si può calcolare il valor medio dell’energia per questa funzione (che in generale non sarà un’autofunzione di H): hHi =

Z

ψ ∗ Hψ dv

(6.36)

dove v rappresenta tutte le coordinate di integrazione. Il principio variazionale afferma che le funzioni ψ per le quali hHi è stazionario—ossia non varia al primo ordine per piccole variazioni di ψ—sono le autofunzioni dell’energia. In altre parole, l’equazione di Schrödinger è equivalente ad una condizione di stazionarietà.

65

6.3.1

Dimostrazione del principio variazionale (I)

Poichè una variazione arbitraria δψ di una funzione d’onda in generale ne distrugge la normalizzazione, è conveniente utilizzare la definizione pi` u generale di valor medio R ∗ ψ Hψ dv (6.37) hHi = R ∗ ψ ψ dv Modificando la ψ in ψ + δψ, il valor medio diventa (ψ ∗ + δψ ∗ )H(ψ + δψ) dv (ψ ∗ + δψ ∗ )(ψ + δψ) dv R ∗ R R ψ Hψ dv + δψ ∗ Hψ dv + ψ ∗ Hδψ dv R R R ψ ∗ ψ dv + δψ ∗ ψ dv + ψ ∗ δψ dv

R

hHi + δhHi = =

R

Z

=

Z

∗

ψ Hψ dv +

Z

∗

δψ Hψ dv +

∗

ψ Hδψ dv ×

R R ∗ δψ ∗ ψ dv ψ δψ dv 1 R R R − 1 − ψ ∗ ψ dv ψ ∗ ψ dv ψ ∗ ψ dv

(6.38)

dove si sono omessi i termini del secondo ordine in δψ, e si è usata l’approssimazione 1/(1+x) ' 1−x valida per piccoli x. Omettendo nuovamente i termini di ordine superiore al primo: R

δhHi =

δψ ∗ Hψ dv R + ψ ∗ ψ dv

R

ψ ∗ Hδψ ∗ dv R − hHi ψ ∗ ψ dv

R ∗ R δψ ∗ ψ dv ψ δψ dv R + R ∗ . (6.39) ψ ∗ ψ dv ψ ψ dv

I due termini nella parentesi tonda sono l’uno il complesso coniugato dell’altro, e lo stesso vale anche per i primi due poichè H è un operatore hermitiano, e soddisfa quindi a Z Z a∗ Hb dv =

b∗ Ha dv

∗

(6.40)

per qualsiasi coppia di funzioni a e b. Pertanto R

δhHi =

δψ ∗ Hψ dv R + c.c. − hHi ψ ∗ ψ dv

R δψ ∗ ψ dv R + c.c. . ∗

ψ ψ dv

(6.41)

Supponiamo ora che ψ sia tale che hHi sia stazionario rispetto a qualsiasi sua variazione. Sar` a allora δhHi = 0, ossia Z

δψ ∗ [H − hHi] ψ dv + c.c. = 0

(6.42)

per una variazione δψ arbitraria, e questo implica che deve essere [H − hHi] ψ = 0

(6.43)

ovvero ψ è una soluzione dell’equazione di Schrödinger: Hψ = Eψ

66

(6.44)

6.3.2

Dimostrazione del principio variazionale (II)

Un altro modo di dimostrare lo stesso principio, utile in seguito, è basato sul metodo dei moltiplicatori di Lagrange. Il metodo afferma che se si vuole rendere stazionario un integrale I0 mantenendo allo stesso tempo costanti altri integrali I1 . . . Ik , si pu` o porre ! δ I0 +

X

λk Ik

=0

(6.45)

k

dove λk sono costanti da determinare. Nel nostro caso avremo Z

ψ ∗ Hψ dv

(6.46)

ψ ∗ ψ dv

(6.47)

δ(I0 + λI1 ) = 0

(6.48)

I0 = Z

I1 = e quindi porremo

con λ da determinare. Procedendo come nella sezione precedente si ha Z

δI0 = Z

δI1 =

δψ ∗ Hψ dv + c.c.

(6.49)

δψ ∗ ψ dv + c.c.

(6.50)

e quindi la condizione da soddisfare è Z

δ(I0 + λI1 ) =

δψ ∗ [H + λ]ψ dv + c.c. = 0

(6.51)

da cui Hψ = −λψ

(6.52)

ossia il moltiplicatore di Lagrange è uguale, a meno del segno, a un autovalore dell’energia. Nuovamente vediamo che gli stati la cui energia media è stazionaria rispetto a qualsiasi variazione della funzione d’onda sono le soluzioni dell’equazione di Schr¨ odinger.

6.3.3

Energia dello stato fondamentale

Siano ψn le autofunzioni di un hamiltoniano H, a cui sono associate energie En : Hψn = En ψn

(6.53)

Supponiamo che lo stato fondamentale corrisponda a n = 0 e abbia quindi energia E0 . Sia ψ una qualunque altra funzione. Dimostriamo che si ha necessariamente R ∗ ψ Hψ dv hHi = R ∗ ≥ E0 (6.54) ψ ψ dv Per dimostrarlo, pensiamo di sviluppare ψ usando la base delle autofunzioni dell’energia. Ci` o è sempre possibile perchè le autofunzioni dell’energia costituiscono un sistema completo e ortonormale, come dimostrato in 4.1. ψ=

X n

67

cn ψn

(6.55)

Sarà allora

P P 2 2 n |cn | (En − E0 ) n |cn | En P = E + hHi = P 0 2 2 n |cn |

n |cn |

(6.56)

Poichè il secondo termine è positivo o nullo, essendo per definizione di stato fondamentale En ≥ E0 , la (6.54) è dimostrata. Questo risultato è semplice ma estremamente importante: ci dice che data una qualsiasi ψ, il suo valor medio dell’energia è sempre una stima superiore dell’energia dello stato fondamentale. Se lo stato fondamentale non è noto, si pu` o quindi pensare di cercare una sua approssimazione facendo variare ψ nell’ambito di un insieme di funzioni di prova e cercando quella funzione che minimizza hHi. Questa è l’essenza del metodo variazionale.

6.3.4

Il metodo variazionale in pratica

Si identifica una famiglia di funzioni d’onda di prova ψ(v; α1 , . . . , αr ), dove v è l’insieme delle variabili, e gli αi sono parametri. L’autovalore dell’energia sarà una funzione del parametri: Z

E(α1 , . . . , αr ) =

ψ ∗ Hψ dv

(6.57)

Il metodo variazionale consiste nel cercare il minimo di E rispetto a variazioni dei parametri, imponendo cioè ∂E ∂E = ... = =0 ∂α1 ∂αr

(6.58)

La ψ che soddisfa a queste condizioni con l’energia pi` u bassa è quella che pi` u si avvicina allo stato fondamentale. Può essere considerata come la miglior approssimazione possibile allo stato fondamentale tra l’insieme delle funzioni di prova. ` chiaro che la scelta della famiglia delle funzioni di prova gioca un ruolo E cruciale e va effettuata con attenzione.

6.4

Problema secolare

Il metodo variazionale pu` o essere ricondotto ad un problema algebrico immaginando di sviluppare la funzione d’onda in una base finita di funzioni, e applicando il metodo variazionale per trovare i coefficienti ottimali dello sviluppo. Basandoci sulla (6.45), ci` o significa calcolare il funzionale (ossia una funzione di funzione) G[ψ] = hψ|H|ψi − hψ|ψi Z

=

ψ ∗ Hψ dv −

Z

ψ ∗ ψ dv

(6.59)

e imporre che G[ψ] sia minimo. Questo dà luogo a una equazione per i coefficienti dello sviluppo che ora determineremo. ` fondamentale notare che la nostra base è costituita da un numero finito E N di funzioni, e quindi non costituirà un sistema completo: ossia non sarà in

68

generale possibile sviluppare una qualsiasi funzione ψ in questa base, tra cui in generale anche le soluzioni esatte dell’equazione di Schrödinger. Quello che faremo è quindi trovare la ψ che meglio si avvicina al vero stato fondamentale nell’ambito di tutte le funzioni esprimibili come combinazione lineare delle N funzioni di base scelte.

6.4.1

Sviluppo in funzioni ortonormali

Supponiamo di avere a disposizione una base di N funzioni bi fra loro ortonormali: Z hbi |bj i ≡ b∗i bj dv = δij (6.60) e sviluppiamo una generica ψ in questa base: ψ=

N X

ci bi

(6.61)

i=1

Sostituendo la (6.61) nella (6.59) si vede immediatamente che quest’ultima prende la forma G(c1 , . . . , cN ) =

X

=

X

c∗i cj Hij −

ij

X

c∗i cj δij

ij

c∗i cj (Hij − δij )

(6.62)

ij

dove si è posto Hij = hbi |H|bj i =

Z

b∗i Hbj dv

(6.63)

Poichè sia H che la base sono dati, Hij è una matrice quadrata di numeri perfettamente nota, e che per la proprietà di hermiticità dell’operatore hamiltoniano tale che Hji = Hij∗ (6.64) (quindi simmetrica nel caso in cui tutti gli elementi siano reali). Come richiesto dal metodo variazionale, minimizziamo la (6.62) rispetto ai coefficienti: ∂G =0 ∂ci

(6.65)

e questo fornisce1 X

(Hij − δij )cj = 0

(6.66)

j 1

Lo scettico pu` o separare i coefficienti in una parte reale e una immaginaria ck = xk + iyk , richiedere che siano nulle sia le derivate rispetto a xk che quelle rispetto a yk , e otterr` a (sfruttando l’hermiticit` a) un sistema

dove Wk =

P j

Wk + Wk∗

=

0

−iWk + iWk∗

=

0

(Hkj − δkj )cj , che ammette come soluzione solo Wk = 0.

69

Notiamo che se la base fosse un sistema completo (e quindi infinita), questa sarebbe una forma dell’equazione di Schrödinger. Abbiamo quindi dimostrato che queste stesse equazioni, nel caso in cui la base sia finita, costituiscono la migliore approssimazione possibile alla vera soluzione secondo il principio variazionale. La (6.66) è un sistema di N equazioni algebriche lineari e omogenee (non ci sono termini costanti) per le N incognite cj . In generale questo sistema ha come unica soluzione possibile tutti i cj nulli (caso che ovviamente non corrisponde ad alcuna funzione d’onda). Per avere soluzioni non nulle è necessario che il determinante dei coefficienti sia nullo: det |Hij − δij | = 0

(6.67)

Ciò corrisponde in pratica a dire che una delle equazioni è una combinazione lineare delle altre, e quindi il sistema si riduce in realtà a un sistema di N − 1 equazioni con N incognite, che ammette soluzione non nulla. La (6.67) è detta equazione secolare. Si tratta di una equazione algebrica di grado N in (come subito si vede sviluppando il determinante e notando che la diagonale principale genera un termine contenente N , e tutte le altre diagonali termini con potenze inferiori), che possiede quindi N radici. Queste radici sono dette gli autovalori. La (6.66) può anche essere scritta in forma matriciale Hc = c

(6.68)

dove H è qui la matrice N × N costituita dagli Hij , e c è un vettore costituito dai ci disposti in colonna. Le soluzioni c sono quindi anche chiamati autovettori. Per ogni radice (autovalore) vi sarà un corrispondente autovettore (determinato a meno di una costante moltiplicativa, fissata dalla normalizzazione). Avremo quindi N autovettori. Si potr` a allora scrivere che vi sono N soluzioni ψk =

X

Cik bi ,

k = 1, . . . , N

(6.69)

i

dove Cik è una matrice costruita disponendo in colonna, fianco a fianco, gli N autovettori, e tali che Hψk = k ψk (6.70) ovvero, in forma matriciale, prendendo la componente i−ima, (Hψk )i =

X

Hij Cjk = k Cik

(6.71)

j

La (6.68) è una equazione comune nell’algebra lineare, ed esistono metodi standard per risolverla. Data una matrice H, si ottengono quindi facilmente— attraverso routine di libreria—la matrice C e un vettore di autovalori . Il processo di risoluzione è anche noto come diagonalizzazione. Questo nome deriva dalla seguente importante proprietà di C. La (6.69) può essere vista come una trasformazione delle N funzioni di partenza in un altro insieme di N funzioni attraverso una matrice di trasformazione. Si può far vedere che

70

se le bi sono fra loro ortonormali anche le ψk lo sono. Si dice allora che la trasformazione è unitaria. Ci` o corrisponde ad affermare che X

∗ Cij Cik = δjk

(6.72)

i

o in notazione matriciale † ∗ (C −1 )ij = Cji ≡ Cij

(6.73)

ossia la matrice inversa è uguale alla trasposta coniugata, ovvero alla matrice aggiunta (cioè C è una matrice unitaria). Consideriamo ora il prodotto di matrici C −1 HC e calcoliamo un suo elemento: (C −1 HC)kn =

X

=

X

(C −1 )ki Hij Cjn

ij ∗ Cik

i

=

X

Hij Cjn

j

X

∗ Cik n Cin

i

= n

X

∗ Cin Cik

i

(6.74)

= n δkn

avendo fatto uso dei risultati precedenti. Si dice allora che la trasformazione C riduce H in forma di una matrice diagonale, i cui N elementi non nulli sono gli autovalori. Possiamo vedere quindi il nostro problema agli autovalori come quello della ricerca di una trasformazione che porti dalla base originale ad una nuova base in cui l’operatore H ha una forma diagonale, ossia agisce sugli elementi della base semplicemente moltiplicandoli per una costante (equazione di Schr¨ odinger).

6.4.2

Sviluppo in funzioni non ortonormali

I metodi di algebra di lineare permettono di trattare agevolmente anche il caso in cui la base è costituita da funzioni non ortonormali tra loro, in cui cioè Sij = hbi |bj i =

Z

b∗i bj dv

(6.75)

non è pari a δij . Gli Sij vengono detti integrali di overlap, per ovvi motivi. Talvolta risulta comodo lavorare con basi di questo tipo. Considerazioni simili a quelle effettuate all’inizio della sezione precedente indicano che in questo caso la (6.62) assume la forma pi` u generale G(c1 , . . . , cN ) =

X

c∗i cj (Hij − Sij )

(6.76)

ij

e la condizione di minimo (6.66) diventa X

(Hij − Sij )cj = 0

j

71

(6.77)

o in forma matriciale Hc = Sc

(6.78)

noto come problema agli autovalori generalizzato. Risolvere un problema agli autovalori generalizzato corrisponde dal punto di vista numerico a risolvere due problemi agli autovalori semplici. Supponiamo infatti di procedere in due stadi, occupandosi prima del problema ausiliario Sd = σd

(6.79)

Questo è del tutto analogo al problema (6.68). Potremo cos`ı trovare una matrice D (ottenuta disponendo gli autovettori per colonne), tale che D−1 SD sia diagonale, e i cui elementi non nulli siano gli autovalori σ. Si otterrà una equazione analoga alla (6.74): X

∗ Dik

X

i

Sij Djn = σn δkn

(6.80)

j

Supponiamo per` o di definire un’altra matrice di trasformazione Dij Aij ≡ √ σj

(6.81)

Si avr` a allora X

A∗ik

i

X

Sij Ajn = δkn

(6.82)

j

o in forma matriciale A∗T SA = I

(6.83)

dove T indica trasposizione. Una matrice A con questa proprietà può dunque essere ottenuta risolvendo un normale problema agli autovalori. Poniamo ora c = Av (6.84) Con questa posizione, l’equazione (6.78) diventa HA v = SA v

(6.85)

Moltiplichiamo a sinistra per A∗T : A∗T HA v = A∗T SA v = v

(6.86)

Pertanto, se risolviamo ora il problema secolare per l’operatore A∗T HA troveremo gli autovalori desiderati per l’energia. Per ottenere gli autovettori nella nostra base di partenza baster` a, secondo la (6.84), applicare l’operatore A su ciascun autovettore.

72

6.5

Programma: hydrogen gauss

Il programma hydrogen gauss.f902 (oppure hydrogen gauss.c3 ) risolve il problema secolare per l’atomo di idrogeno utilizzando due diverse basi non ortonormali: 1. una base gaussiana “onda S”: 2

bi (r) = e−αi r ;

(6.87)

2. una base gaussiana “onda P”, di cui esistono tre scelte possibili, corrispondenti a valori diversi di m: 2

2

bi (r) = xe−αi r ,

bi (r) = ye−αi r ,

2

bi (r) = ze−αi r .

(6.88)

L’operatore hamiltoniano del problema è ovviamente H=−

¯h2 ∇2 Zqe2 − 2me r

(6.89)

e per l’atomo di idrogeno Z = 1. I calcoli per l’onda S e l’onda P sono del tutto indipendenti: si tratta di due basi distinte. La base in onda P è chiaramente inadatta a descrivere lo stato fondamentale, perchè non ha la simmetria (angolare) giusta, e viene inclusa a scopo didattico. Il codice legge da file una lista di esponenti αi e procede quindi a valutare tutti gli elementi delle matrici Hij e Sij . Il calcolo è basato sulle espressioni degli integrali calcolati analiticamente. In particolare si trova per l’onda S Z

Sij =

−(αi +αj )r2 3

e

d r=

π αi + αj

!3/2

(6.90)

e inoltre i termini cinetico e coulombiano di Hij sono rispettivamente HijK

Z

=

−αi r2

e

HijV

"

#

¯h2 ∇2 −αj r2 3 ¯h2 6αi αj − e d r= 2me 2me αi + αj Z

=

−αi r2

e

"

π αi + αj

!3/2

(6.91)

#

Zq 2 2πZqe2 2 − e e−αj r d3 r = − r αi + αj

(6.92)

Per l’onda P si procede analogamente, utilizzando le corrispondenti espressioni per gli integrali. Il codice procede dunque a chiamare la subroutine diag che risolve il problema secolare generalizzato (ossia applica il principio variazionale), ritornando il vettore e contenente gli autovalori in ordine crescente di energia, e la matrice v contenente gli autovettori, ossia i coefficienti dello sviluppo. 2 3

http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/F90/hydrogen gauss.f90 http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/C/hydrogen gauss.c

73

Al suo interno, diag effettua il calcolo in due stadi descritto nella sezione precedente, delegando la risoluzione del problema agli autovalori semplice alla subroutine dsyev della libreria di algebra lineare LAPACK4 , ed utilizzando alcune chiamate addizionali a subroutines della libreria BLAS5 . Informazioni addizionali su queste librerie sono disponibili a partire dalla pagina del corso lapack.html6 . Notare come dopo la prima diagonalizzazione della matrices S, diag getta via tutti gli autovettori corrispondenti ad autovalori molto piccoli, entro una soglia di tolleranza numerica. Questi corrispondono a combinazioni quasi linearmente dipendenti dalle altre, che quindi implicano una quasi singolarità della matrice. Il numero di vettori di base linearmente indipendenti trovati viene riportato nell’output. Tale procedura non è in linea di principio necessaria: esistono routines di LAPACK che risolvono il problema secolare generalizzato Hψ = Sψ in una sola chiamata. Se la base è ben scelta, la matrice S non ha autovalori molto piccoli e le conseguenti instabilità numeriche. Tuttavia è facile verificare che in presenza di molte funzioni di base, o di funzioni di base mal scelte (per esempio troppo vicine) si possano avere autovalori quasi singolari. Una buona scelta dei coefficienti αj è fondamentale per avere un’alta accuratezza e assenza di instabilit` a. Il programma procede quindi a scrivere i coefficienti dello sviluppo nel file s-coeff.out (p-coeff.out): per ogni funzione j, e la funzione d’onda dello stato fondamentale nel file s-wfc.out (p-wfc.out).

6.5.1

Laboratorio

• osservare e discutere lo stato fondamentale ottenuto con la base d’onde P • verificare l’accuratezza degli autovalori dell’energia • osservare gli effetti legati al numero di funzioni di base, e alla scelta dei parametri α. Provate per esempio a scegliere le lunghezze caratteristiche √ delle gaussiane, λ = 1/ α, uniformemente distribuite fra un λmin e un λmax opportunamente scelto. • confrontare la soluzione con quella ottenuta dal programma hydrogen. Provare in particolare la seguente base ottimizzata di quattro gaussiane: α1 = 0.121949, α2 = 0.444529, α3 = 1.962079, α4 = 13.00773 (a.u.). • Per Z > 1, come riscalereste i coefficienti delle gaussiane adatti per Z = 1?

4

http://www.netlib.org/lapack/ http://www.netlib.org/blas/ 6 http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/lapack.html 5

74

6.6

Base di onde piane

Un altro tipo di base comunemente impiegato in fisica è la base di onde piane. Questa è strettamente legata al concetto di trasformata di Fourier. Una funzione f (x) definita su tutto l’asse reale può essere sviluppata in componenti di Fourier, fe(k): ∞ 1 √ fe(k)eikx dk 2π −∞ Z ∞ 1 √ f (x)e−ikx dx. 2π −∞

Z

f (x) = fe(k) =

(6.93) (6.94)

Per un intervallo finito [−L/2, L/2], possiamo invece scrivere f (x) = fe(kn ) =

1 Xe √ f (kn )eikn x L n 1 √ L

Z L/2

(6.95)

f (x)e−ikn x dx

(6.96)

−L/2

con kn = 2πn/L, n = 0, ±1, ±2, .... Da notare che la f (x) dell’Eq.6.95 è per costruzione una funzione periodica di periodo L: vale f (x + L) = f (x), come si verifica immediatamente. Ci` o implica che deve valere f (−L/2) = f (+L/2) (condizioni al contorno periodiche). Le formule qui sopra si generalizzano in modo banale a tre o pi` u dimensioni. Nel seguito ci limiteremo ad un caso unidimensionale. Definiamo quindi la nostra base di onde piane bi (x) secondo la (6.95): 1 bi (x) = √ eiki x , L

ki =

2π i, L

i = 0, ±1, ±2, ..., ±N

(6.97)

e i relativi coefficienti ci della funzione d’onda ψ(x) come Z L/2

ci =

−L/2

b∗i (x)ψ(x)dx = hbi |ψi,

ψ(x) =

X

ci bi (x).

(6.98)

i

Questa base, composta da 2N + 1 funzioni, diventa completa nel limite N → ∞ (è una conseguenza delle note proprietà delle serie di Fourier). E’ immediato verificare che questa base è ortonormale: Sij = hbi |bj i = δij . La soluzione del problema di una particella sotto ad un potenziale richiede quindi la diagonalizzazione della matrice hamiltoniana, i cui elementi di matrice: Hij = hbi |H|bj i = hbi |

p2 + V (x)|bj i 2m

(6.99)

sono banali. Il termine cinetico è diagonale: hbi |

p2 ¯h2 |bj i = − 2m 2m

Z L/2 −L/2

b∗i (x)

d 2 bj ¯h2 ki2 (x)dx = δ . ij dx2 2m

(6.100)

Il termine di potenziale non è altro che una trasformata di Fourier: 1 hbi |V (x)|bj i = L

Z L/2

1 V (x)e−i(ki −kj )x dx = √ Ve (ki − kj ). L −L/2

75

(6.101)

Di conseguenza, gli elementi di matrice del potenziale tenderanno a zero per grandi valori di ki − kj tanto pi` u rapidamente quanto pi` u il potenziale varia lentamente, e viceversa. Potenziali (e funzioni d’onda) che variano su di una scala tipica di lunghezza λ avranno trasformata di Fourier significativa fino a kmax ∼ 2π/λ. In questo modo possiamo stimare il numero di onde piane necessarie per risolvere un problema.

6.7

Programma: pwell

Consideriamo un semplice problema: i livelli di una buca di potenziale a a per x < − , x > 2 2 a a V (x) = −V0 per − ≤ x ≤ 2 2 V (x) = 0

(6.102) (6.103)

dove V0 > 0, a < L. Gli elementi di matrice dell’hamiltoniana sono dati dal’Eq. (6.100) per la parte cinetica, dall’Eq.(6.101) per la parte potenziale. Quest’ultima si calcola esplicitamente: 1 hbi |V (x)|bj i = − L

Z a/2 −a/2

V0 e−i(ki −kj )x dx

(6.104)

a/2

V0 e−i(ki −kj )x = − V0 L −i(ki − kj ) −a/2 =

V0 sin (a(ki − kj )/2) , L (ki − kj )/2

(6.105) ki 6= kj .

(6.106)

Il caso ki = kj va trattato a parte e dà Ve (0) =

V0 a . L

(6.107)

Il programma pwell.f90 genera i ki , riempe la matrice Hij e la diagonalizza. Il programma usa unit` a nelle quali h ¯ 2 /2m = 1. I dati di input sono: larghezza (a) e profondit` a (V0 ) della buca, larghezza della scatola (L), numero di onde piane (2N + 1). In uscita, il programma stampa i tre livelli pi` u bassi e la funzione d’onda dello stato fondamentale (su file).

6.7.1

Laboratorio

• Osservate la convergenza dei risultati con il numero di onde piane e verificate la forma della funzione d’onda. Verificate l’energia (per esempio: per V0 = 1, a = 2, il risultato esatto è E = 0.4538) in particolare nel caso limite V0 → ∞. • Osservate la convergenza rispetto a L. Si noti come per valori di L non molto grandi, l’energia calcolata con il metodo variazionale è inferiore al valore esatto. Perché?

76

• Provate a modificare il programma per una buca di forma gaussiana (la cui trasformata di Fourier è nota in forma analitica). A parità di ”larghezza”, converge pi` u rapidamente il problema a buca quadrato o quello a buca gaussiana? • Sappiamo che per un potenziale simmetrico rispetto alla parità, V (−x) = V (x), le soluzioni hanno parità alternata: pari lo stato fondamentale, dispari il primo stato eccitato, e cos`ı via. Sfruttate questo fatto per ridurre il problema in due sottoproblemi, uno per gli stati pari e uno per gli stati dispari, usando funzioni seno e coseno ottenute come combinazioni delle onde piane usate finora (attenzione alla normalizzazione corretta e al termine kn = 0). Cosa ci si guadagna in termini computazionali?

77

Chapter 7

Atomi a pi` u elettroni 7.1

Lo spin

Gli esperimenti indicano che alle particelle si deve associare un “momento angolare intrinseco”, o spin, indipendentemente dalla loro natura (particelle elementari o con una struttura interna), e indipendentemente dal loro moto nello spazio. Si tratta di una variabile in pi` u necessaria per descrivere il sistema. Lo spin di una particella pu` o comparire nell’hamiltoniana sia perchè accoppiato ad esempio con un campo magnetico esterno, o perchè accoppiato con lo stesso momento angolare orbitale, e quindi dare origine a effetti chiaramente osservabili. Lo spin S è un momento angolare—nel senso che soddisfa all’algebra di commutazione (A.39)—e quindi valgono tutti i risultati della sez.A.7. Inoltre, l’autovalore di S 2 è una proprietà intrinseca di un tipo di particelle, e quindi ciò che varia è la proiezione ms dello spin lungo un asse. In particolare lo spin dell’elettrone è 1/2. Ciò significa che S 2 ha il valore fissato (3/4)¯ h2 , e la sua proiezione lungo un asse ha autovalori −¯h/2 o +¯ h/2. Lo spin di un elettrone pu` o essere rappresentato dalle matrici di Pauli: 



σ ¯h  x  S =  σ y  , σx = 2 σz

0 1 1 0

!

, σy =

0 −i i 0

!

, σz =

1 0 0 −1

!

. (7.1)

Tali matrici agiscono su vettori di dimensione 2, gli spinori: φ=

α β

!

(7.2)

,

in cui la parte superiore (α) rappresenta lo stato con spin up, la parte inferiore (β) lo stato con spin down. Si verifica facilmente che le matrici di Pauli danno le corrette relazioni di commutazione, e in pi` u 3¯h2 ¯h2 2 S = (σx + σy2 + σz2 ) = 4 4 2

78

1 0 0 1

!

=

3¯h2 . 4

(7.3)

7.2

Composizione di momenti angolari: la rappresentazione accoppiata

Consideriamo un sistema in cui J1 e J2 sono due operatori momento angolare che commutano tra loro. Ciò accade quando si riferiscono a sistemi fisici indipendenti; ad esempio i momenti angolari di due particelle diverse, oppure il momento angolare orbitale e quello di spin di una stessa particella nell’assunzione che non vi siano interazioni che li accoppiano. Avremo quindi quattro osservabili commutanti per descrivere il sistema: J12 , J1z , J22 e J2z , e gli autostati comuni saranno caratterizzati dall’insieme di numeri quantici j1 , m1 , j2 e m2 . Dati j1 e j2 , avremno quindi (2j1 + 1)(2j2 + 1) stati distinti. Esiste anche un’altro utile insieme di osservabili per descrivere lo stesso sistema. Definiamo l’operatore momento angolare totale J = J1 + J2

(7.4)

` immediato verificare che anche J deve soddisfare all’algebra di commutazione E del momento angolare. Infatti [Jx , Jy ] = = = =

[J1x + J2x , J1y + J2y ] [J1x , J1y ] + [J2x , J2y ] i¯hJ1z + i¯hJ2z i¯hJz

(7.5)

avendo sfruttato il fatto che i commutatori fra componenti relative al sistema 1 e al sistema 2 sono per ipotesi nulle. Avremo quindi [Jz , J 2 ] = 0. Possiamo descrivere allora il sistema anche usando i quattro operatori J12 , 2 J2 , J 2 e Jz . Si tratta della cosiddetta rappresentazione accoppiata (il motivo di questo nome sar` a chiaro tra breve). Per dimostrare che commutano tutti fra loro ci resta da vedere che [J12 , J 2 ] = [J22 , J 2 ] = 0, e che [J12 , Jz ] = [J22 , Jz ] = 0. [J12 , J 2 ] = = = = = e

[J12 , Jx Jx + Jy Jy + Jz Jz ] Jx [J12 , Jx ] + [J12 , Jx ]Jx + (. . .) Jx [J12 , J1x + J2x ] + [J12 , J1x + J2x ]Jx + (. . .) Jx [J12 , J1x ] + [J12 , J1x ]Jx + (. . .) 0 [J12 , Jz ] = [J12 , J1z + J2z ] = [J12 , J1z ] = 0

(7.6)

(7.7)

Indichiamo con j1 , j2 , j e m i numeri quantici che caratterizzano gli autovalori dei nostri operatori. Dati j1 e j2 , j varierà da un certo jmin a un certo jmax , che ora identificheremo. m è la proiezione di Jz = J1z + J2z , e quindi per definizione dovrà essere m = m1 + m2 . Da questo si desume che il massimo valore possibile per m è j1 + j2 , ma il massimo valore possibile per m è anche pari a jmax . Dunque jmax = j1 + j2 .

79

jmin si ottiene imponendo che il numero totale di stati sia lo stesso ottenuto nella prima rappresentazione: j1X +j2

(2j + 1) = (2j1 + 1)(2j2 + 1)

(7.8)

j=jmin

Si pu` o verificare che questo fornisce jmin = |j1 − j2 |. Dunque j = |j1 − j2 |, . . . , j1 + j2 . Si può pensare che nel caso j = |j1 − j2 | i due vettori abbiano la stessa direzione ma verso opposto, e nel caso j = j1 + j2 la stessa direzione e verso; i casi intermedi corrispondono a momenti angolari che puntano in direzioni diverse.

7.2.1

Esempio: singoletti e tripletti

Supponiamo che sia j1 = 1/2 e j2 = 1/2. Avremo quattro stati indipendenti, con m1 = ±1/2 e m2 = ±1/2. Passiamo alla rappresentazione accoppiata. I possibili valori di j sono j = 0 e j = 1. Nel caso j = 0 si può avere solo m = 0 (“stato di singoletto”). Nel caso j = 1 si avrà m = 0, ±1 (“stati di tripletto”). In totale vi sono sempre quattro stati.

7.2.2

Presenza di accoppiamento

Perchè abbiamo introdotto questa rappresentazione, apparentemente equivalente a quella che considera i due momenti angolari separatamente? Il motivo è che in molti casi sono presenti termini nell’hamiltoniana che accoppiano i momenti angolari tra loro. Comune è ad esempio il caso di una “coppia” che tende ad allineare i due vettori: H = . . . − AJ1 · J2 (7.9) In presenza di tale termine, J1z e J2z non sono pi` u conservati, ossia non sono pi` u buoni numeri quantici. Infatti [J1z , −AJ1 · J2 ] = −A[J1z , J1x J2x + J1y J2y + J1z J2z ]

(7.10)

= −A[J1z , J1x ]J2x − A[J1z , J1y ]J2y

(7.11)

= −i¯hA(J1y J2x − J1x J2y )

(7.12)

e questo operatore in generale non è nullo. Viceversa, è immediato vedere che J1z +J2z è conservato, in quanto [J2z , −AJ1 · J2 ] risulta essere uguale al termine calcolato qua sopra con segno opposto.

7.3

Particelle identiche: principio di indistinguibilit` a

Nella meccanica quantistica non esiste il concetto di traiettoria della meccanica classica, che presuppone la conoscenza simultanea della posizione e della velocità delle particelle. Ci` o ha delle importanti implicazioni. Supponiamo di considerare due particelle del tutto identiche (ad esempio due elettroni), e di determinare con elevata precisione la loro posizione ad un

80

certo istante t, trovando due posizioni r1 e r2 . Supponiamo di ripetere la misura ad un successivo istante t0 , trovando delle posizioni r01 e r02 . Siamo in grado di dire se la particella in r01 era quella che si trovava in r1 , oppure quella che si trovava in r2 all’osservazione precedente? La risposta è no. Questo è un principio generale che prende il nome di “principio di indistinguibilit` a”: Dato un sistema contenente N particelle fra loro identiche, è impossibile che una misura dia risultati diversi se si immagina di scambiare fra loro due particelle. In altre parole, il sistema deve essere simmetrico rispetto a tutte le permutazioni possibili. Immaginiamo per il momento di avere a che fare con un sistema in cui le particelle non interagiscono fra loro. Il problema di Schrödinger è allora separabile in N equazioni ad una particella, ed è possibile scrivere una soluzione per la funzione d’onda complessiva in forma di prodotto di soluzioni delle funzioni d’onda per le singole particelle. Per un sistema a due particelle: ψ(1, 2) = φ1 (1)φ2 (2)

(7.13)

dove abbiamo indicato con (1) e (2) le variabili associate alle due particelle (tipicamente per ogni particella le tre coordinate di posizione e la variabile di spin intrinseco), mentre invece gli indici bassi in φ1 e φ2 indicano la particolare funzione d’onda scelta, classificata coi numeri quantici del problema a una particella. Ebbene, la (7.13) non è una soluzione accettabile perchè—pur soddisfacendo all’equazione di Schr¨ odinger—viola il principio di instinguibilità. Lo scambio delle particelle porta infatti ad una funzione ψ(1, 2) = φ2 (1)φ1 (2)

(7.14)

che pure è soluzione dell’equazione di Schrödinger ed è nettamente diversa dalla precedente. Ad esempio φ1 potrebbe essere un orbitale 1s, φ2 un orbitale 2p, e le due funzioni sopra darebbero origine a distribuzioni di densità (quantità misurabili) diverse. ` comunque possibile costruire soluzioni che soddisfano al principio di inE distinguibilit` a combinando opportunamente le (7.13) e (7.14): 1 ψs (1, 2) = √ [φ1 (1)φ2 (2) + φ2 (1)φ1 (2)] 2

(7.15)

(funzione d’onda simmetrica) e 1 (7.16) ψa (1, 2) = √ [φ1 (1)φ2 (2) − φ2 (1)φ1 (2)] 2 √ (funzione d’onda antisimmetrica) I fattori 1/ 2 servono a mantenere le normalizzazioni corrette. Il principio di indistinguibilità è ovviamente soddisfatto per la (7.15). Nel caso della (7.16), lo scambio delle particelle porta ad un cambiamento di segno, ma le quantità osservabili (associate a |ψ|2 ) restano inalterate. Quale delle due trasformazioni va scelta? La risposta dipende dal tipo di particella, come discusso nella prossima sezione.

81

7.4

Operatori di permutazione

Consideriamo un sistema a due particelle e chiamiamo P l’operatore che le scambia. Ossia: (7.17) P ψ(1, 2) ≡ ψ(2, 1) per qualsiasi funzione d’onda. Si vede subito che P è un ”operatore idempotente”, ossia soddisfa a P 2 = 1, e quindi P = P −1 . Il principio di indistinguibilità ci dice che qualsiasi misura deve dare lo stesso risultato se effettuata sullo stato ψ o sullo stato P ψ. Se l’operatore A corrisponde ad una generica osservabile e il sistema è in un autostato di A: Aψ = aψ

(7.18)

(ossia, una misura di A d` a come risultato un numero a ben definito), allora deve anche essere vero che (7.19) AP ψ = aP ψ D’altra parte applicando l’operatore P a sinistra e a destra nella prima equazione si ha anche (7.20) P Aψ = aP ψ Sottraendo tra loro le ultime due equazioni e notando che questo deve valere per qualsiasi stato fisico ψ, deve allora essere [P, A] = 0

(7.21)

ossia l’operatore P deve commutare con qualsiasi osservabile fisica, inclusa l’energia: (7.22) [P, H] = 0 e quindi è una quantit` a conservata. Dalla proprietà di idempotenza si inferisce inoltre che se P ψ = λψ, deve valere λ2 = 1, ovvero i suoi autovalori possono essere solamente +1 oppure −1. Risulta che il segno dell’autovalore dell’operatore di scambio, o parit` a, è una propriet` a intrinseca del tipo di particella. Le particelle si dividono in • bosoni: P ψ = +ψ • fermioni: P ψ = −ψ dove P è riferito a una qualsiasi coppia di particelle di quel tipo. Un insieme di particelle si comporta dunque sempre in un dato modo, che dipende esclusivamente dal carattere di bosone o fermione della particella.

7.5

Caso di pi` u particelle e sistemi composti

Se un sistema è costituito da N particelle identiche anzichè da due, i risultati ottenuti si generalizzano facilmente. Indicato con P un generico operatore di

82

permutazione che applicato su uno stato di N particelle dà lo stato equivalente in cui le particelle sono state tra loro permutate in modo arbitrario, si avrà Pψ = ψ

(7.23)

Pψ = (−1)M ψ

(7.24)

per qualunque sistema di bosoni, e

per un sistema di fermioni, dove M è il numero di scambi di coppie necessario per arrivare dallo stato iniziale a quello finale. Da ci` o segue subito anche la regola per trovare il carattere di un sistema costituito da particelle non elementari, ma internamente costituite da k bosoni e ` fermioni. Poichè scambiare tra loro due di queste particelle significa scambiare tra loro k bosoni e ` fermioni, si ha subito che la particella composta è un bosone se ` è pari, e un fermione se ` è dispari, indipendentemente da k. Con argomenti di meccanica quantistica relativistica si può dimostrare che • le particelle a spin intero sono bosoni • le particelle a spin semiintero sono fermioni. ` facile vedere che nel caso di particelle composte la regola di composizione dei E momenti angolari d` a un risultato consistente: l’insieme di k particelle a spin intero e ` particelle a spin semiintero dà una particella a spin intero se ` è pari, o una a spin semiintero se ` e dispari, indipendentemente da k.

7.6

Determinanti di Slater

Il caso che pi` u ci interessa è quello degli elettroni, che avendo spin 1/2 sono fermioni, e quindi la loro funzione d’onda deve essere antisimmetrica rispetto allo scambio di qualsiasi coppia. Supponiamo che uno stato sia descrivibile in termini di un prodotto di funzioni d’onda ad un elettrone. Come già discusso per due particelle, un semplice prodotto ψ(1, 2, . . . , N ) = φ1 (1)φ2 (2) . . . φN (N ) (7.25) non soddisfa al principio di indistinguibilità perchè non è un autostato degli operatori di permutazione. ` possibile per` E o costruire una soluzione antisimmetrica per scambio in forma di un determinante: φ1 (1) φ1 (2) . . . φ1 (N ) φ2 (1) φ2 (2) . . . φ2 (N ) 1 . . . . ψ(1, 2, . . . , N ) = √ . . . . N ! . . . . φ (1) φ (2) . . . φ (N ) N N N

(7.26)

Scambiare fra loro due particelle equivale infatti a scambiare fra loro due colonne, e per le propriet` a del determinante questo porta ad un cambiamento di segno.

83

Notiamo per` o anche che se due qualsiasi delle righe sono fra loro uguali, il determinante si annulla e quindi una tale funzione d’onda non corrisponde ad alcun stato fisico. Pertanto tutte le φi devono essere diverse; due (o pi` u) fermioni identici non possono trovarsi nello stesso stato. Si tratta del noto principio di esclusione di Pauli.

7.7

Atomi a due elettroni

Supponiamo che lo spin sia separabile dalle coordinate (cosa senz’altro vera se l’hamiltoniano non contiene termini esplicitamente dipendenti dallo spin). In tal caso si potr` a scrivere ψ(1, 2) = Φ(r1 , r2 )χ(σ1 , σ2 )

(7.27)

dove Φ è funzione solo delle coordinate r e χ solo degli spin σ. La ψ(1, 2) è sempre antisimmetrica perchè gli elettroni sono fermioni. Tuttavia, è chiaro che è possibile ottenere questo risultato con una Φ antisimmetrica e una χ simmetrica, oppure con una Φ simmetrica e una χ antisimmetrica. Dati le autofunzioni di spin del singolo elettrone, ciascuna delle quali ha due valori possibili che indichiamo semplicemente con v+ e v− , possiamo costruire tre funzioni simmetriche dello spin: χ1,1 = v+ (σ1 )v+ (σ2 ) 1 χ1,0 = √ [v+ (σ1 )v− (σ2 ) + v− (σ1 )v+ (σ2 )] 2 χ1,−1 = v− (σ1 )v− (σ2 )

(7.28) (7.29) (7.30)

e una antisimmetrica: 1 χ0,0 = √ [v+ (σ1 )v− (σ2 ) − v− (σ1 )v+ (σ2 )] 2

(7.31)

Quelle simmetriche costituiscono un “tripletto” e corrispondono a uno stato del sistema a due elettroni con spin complessivo pari a 1, e tre possibili valori per la sua proiezione lungo z: -1, 0 e +1. Quella antisimmetrica costituisce un “singoletto” e corrisponde a uno stato con spin complessivo 0. Il valore dello spin complessivo determina quindi la simmetria della parte di spin, e di conseguenza quella della parte configurazionale. La funzione d’onda configurazionale antisimmetrica tende a “respingere” i due elettroni, in quanto non permette che essi possano essere vicini (la funzione d’onda tende ad annullarsi quando gli elettroni vengono portati nella stessa posizione). Per effetto della repulsione elettrostatica, ciò fa s`ı che l’energia risultante sia pi` u bassa di quella del corrispondente caso simmetrico, in cui gli elettroni hanno elevata probabilit` a di trovarsi vicini. Per questo motivo, fra gli stati eccitati dell’elio in cui uno dei due elettroni si trova in un orbitale 2s, lo stato in cui i due spin sono allineati (ortoelio, tripletto, parte di spin simmetrica e parte configurazionale antisimmetrica) ha energia pi` u bassa di quello in cui i due spin sono opposti (paraelio, singoletto, parte di spin antisimmetrica e parte configurazionale simmetrica).

84

Il problema non si pone invece per lo stato fondamentale, in cui entrambi gli elettroni si trovano in orbitali 1s e quindi, come discusso nella sezione 7.8, la funzione d’onda configurazionale deve essere simmetrica. Il concetto di “orbitale” verr` a chiarito meglio nel capitolo seguente. E’ lecito chiedersi quale sia la relazione fra la forma della funzione d’onda data dell’Eq.(7.27) e quella ottenuta come determinante di Slater, Eq.(7.26): 1 ψ(1, 2) = √ [φ1 (1)φ2 (2) − φ2 (1)φ1 (2)] . 2

(7.32)

Si pu` o verificare che una funzioni di tripletto con spin massimo può essere scritta come determinante di Slater delle due funzioni φ1 (r)v+ (σ) e φ2 (r)v+ (σ). Si ottiene 1 (7.33) Φa (r1 , r2 ) = √ [φ1 (r1 )φ2 (r2 ) − φ2 (r1 )φ1 (r2 )] . 2 La funzione di tripletto con spin minimo si ottiene in modo analogo, sostituendo v− a v+ nella parte di spin. La funzione di tripletto con spin zero non è invece ottenibili direttamente da un singolo determinante di Slater: bisogna sovrapporne due, uno generato da φ1 (r)v+ (σ) e φ2 (r)v− (σ), l’altro con le parti di spin scambiate. La sovrapposizione ortogonale degli stessi determinanti di Slater dà la funzione di singoletto, con 1 Φs (r1 , r2 ) = √ [φ1 (r1 )φ2 (r2 ) + φ2 (r1 )φ1 (r2 )] . 2

(7.34)

Ovviamente se φ1 = φ2 esiste un solo determinante di Slater che coincide con lo stato di singoletto, e Φs (r1 , r2 ) = φ(r1 )φ(r2 ).

7.8

(7.35)

Trattamento perturbativo dell’atomo di elio

L’atomo di elio è caratterizzato da un operatore hamiltoniano H=−

¯h2 ∇21 Zqe2 ¯h2 ∇22 Zqe2 q2 − − − + e 2me r1 2me r2 r12

(7.36)

dove r12 = |r2 − r1 | è la distanza tra i due elettroni. L’ultimo termine, corrispondente alla repulsione coulombiana tra i due elettroni, li accoppia tra loro e rende il problema non separabile. In prima approssimazione è possibile però considerare l’interazione tra elettroni q2 (7.37) V = e r12 come una perturbazione al problema descritto da H0 = −

¯h2 ∇21 Zqe2 ¯h2 ∇22 Zqe2 − − − 2me r1 2me r2

85

(7.38)

che è facile da risolvere in quanto si separa in due problemi di un singolo elettrone in un campo di forze centrale coulombiano, ossia il problema di un atomo idrogenoide con nucleo di carica Z. Lo stato fondamentale di un singolo elettrone è dato da una funzione del tipo (5.71), che scriviamo in a.u.: Z 3/2 φ0 (ri ) = √ e−Zri π

(7.39)

(”orbitale 1s”). Notiamo che possiamo assegnare ad entrambi gli elettroni la stessa funzione d’onda, purchè il loro spin sia opposto (se il loro spin fosse uguale, uno dei due dovrebbe essere portato in uno stato eccitato altrimenti il principio di esclusione verrebbe violato). La funzione d’onda totale imperturbata è allora semplicemente il prodotto ψ 0 (r1 , r2 ) =

Z 3 −Z(r1 +r2 ) e π

(7.40)

che è gi` a una funzione simmetrica. L’energia del corrispondente stato fondamentale sar` a la somma delle energie dei due atomi idrogenoidi: E0 = −2Z 2 Ry = −8Ry

(7.41)

essendo Z = 2. La repulsione tra elettroni dovrà alzare l’energia, rendendola cioè meno negativa. Nella teoria delle perturbazioni al primo ordine, E − E0 = hψ0 |V |ψ0 i Z Z6 2 −2Z(r1 +r2 ) 3 3 e d r1 d r2 = π2 r12 5 = ZRy 4

(7.42) (7.43) (7.44)

come si ottiene calcolando l’integrale. Per Z = 2 la correzione è pari a 2.5 Ry, e fornisce un’energia E = −8 + 2.5 = −5.5 Ry. Il valore sperimentale è pari a −5.8074 Ry. L’approssimazione perturbativa non è precisa, ma fornisce una stima ragionevole della correzione pur essendo la “perturbazione” in questo caso di notevole entit` a.

7.9

Trattamento variazionale dell’atomo di elio

Un esempio di applicazione del metodo variazionale può essere fornito nuovamente dall’atomo di elio. L’hamiltoniano, già incontrato nella sezione 7.8, è H=−

q2 ¯h2 ∇21 Zqe2 ¯h2 ∇22 Zqe2 − − − + e 2me r1 2me r2 r12

(7.45)

con Z = 2. Se non vi fosse il termine di repulsione fra i due elettroni il problema sarebbe separabile. Per ciascun elettrone si avrebbe uno stato fondamentale (in a.u.) Z 3/2 (7.46) φ(ri ) = √ e−Zri π

86

con energia associata −Z 2 , e il loro prodotto sarebbe lo stato fondamentale del sistema a due elettroni Z 3 −Z(r1 +r2 ) e π

ψ(r1 , r2 ) =

(7.47)

con energia associata −2Z 2 Ry, ossia −8 Ry. Questo scenario implica l’aver assegnato spin opposti ai due elettroni, averli quindi collocati entrambi in uno stato 1s (n = 1, ` = 0), e aver costruito una funziona d’onda complessiva antisimmetrica in cui la parte di spin è antisimmetrica e quella orbitale simmetrica. L’effetto di ciascun elettrone sull’altro sarà quello di schermare parzialmente il nucleo. Per tener conto della repulsione fra elettroni, possiamo pensare di adottare delle funzioni di prova del tipo (7.47), in cui però sostituiamo la vera carica del nucleo Z con una “carica efficace” Ze , che ci aspettiamo essere pi` u piccola di Z. Questo sar` a il parametro che cercheremo di ottimizzare in modo variazionale. Assumiamo dunque ψ(r1 , r2 ; Ze ) =

Ze3 −Ze (r1 +r2 ) e π

(7.48)

e riscriviamo cos`ı l’hamiltoniano: "

#

"

#

q2 ¯h2 ∇21 Zqe2 ¯h2 ∇22 Zqe2 (Z − Ze )qe2 (Z − Ze )qe2 − − − − + e H= − + − 2me r1 2me r2 r1 r2 r12 (7.49) Calcoliamo ora Z

E(Ze ) =

ψ ∗ (r1 , r2 ; Ze )Hψ(r1 , r2 ; Ze ) d3 r1 d3 r2

(7.50)

Il contributo all’energia dovuto alla prima parentesi quadra in (7.49) è −2Ze2 : si tratta infatti di un problema idrogenoide con nucleo di carica Ze e due elettroni non interagenti. Sviluppando gli integrali rimanenti (e notando che due di essi sono uguali per simmetria) sar` a quindi (in a.u.) E(Ze ) = −2Ze2 −

Z

|ψ|2

4(Z − Ze ) 3 3 d r1 d r2 + r1

Z

|ψ|2

2 3 3 d r1 d r2 r12

(7.51)

con

Ze6 −2Ze (r1 +r2 ) e π2 Gli integrali possono essere calcolati e il risultato è |ψ|2 =

5 27 E(Ze ) = −2Ze2 − 4(Z − Ze )Ze + 2 Ze = 2Ze2 − Ze 8 4

(7.52)

(7.53)

(dove si è esplicitamente inserito Z = 2). La minimizzazione di E(Ze ) rispetto a Ze porta immediatamente a 27 = 1.6875 16

(7.54)

729 = −5.695 Ry 128

(7.55)

Ze = e corrispondentemente E=−

87

Il risultato è decisamente migliore di quello (−5.50 Ry) ottenuto col metodo perturbativo, anche se esiste ancora una discrepanza non trascurabile rispetto al valore sperimentale (−5.8074 Ry). Naturalmente è possibile migliorare il risultato variazionale adottando funzioni di prova pi` u ricche. Questo è quanto viene effettuato dal metodo di Hartree-Fock descritto nel prossimo capitolo, in cui la funzione d’onda complessiva è ancora scritta come un prodotto di funzioni a un elettrone, le quali vengono per` o ottimizzate (ossia non sono semplici esponenziali). Risultati ancora migliori possono essere ottenuti mediante funzioni di prova pi` u complesse di un semplice prodotto. Per esempio, supponiamo di cercare una funzione d’onda radiale del tipo ψ(r1 , r2 ) = [f (r1 )g(r2 ) + g(r1 )f (r2 )] ,

(7.56)

dove le due funzioni f e g sono funzioni idrogenoidi come nella 7.46, ma possono avere valori di Z diversi. Minimizzando rispetto ai due parametri Zf e Zg , si trova Zf = 2.183, Zg = 1.188, e un’energia di E = −5.751 Ry, cioè una riduzione di oltre il 50% dell’errore rispetto al caso di una sola Z effettiva. Si noti come le due funzioni non sono affatto simili!

7.10

Programma: helium gauss

Il programma helium gauss.f901 (oppure helium gauss.c2 ) ricerca lo stato fondamentale dell’atomo di elio, usando lo sviluppo in base di funzioni gaussiane già introdotte per il programma hydrogen gauss. Si cerca la soluzione come prodotto di una parte orbitale simmetrica e di una parte di spin antisimmetrica avente spin totale S = 0. La parte orbitale è sviluppata in funzioni prodotto simmetrizzate di gaussiane, Bk : ψ(r1 , r2 ) =

X

ck Bk (r1 , r2 ).

(7.57)

k

Nel caso in cui le gaussiane, bi , sono di tipo S, come in Eq.(6.87), avremo: 1 Bk (r1 , r2 ) = √ bi(k) (r1 )bj(k) (r2 ) + bi(k) (r2 )bj(k) (r1 ) 2

(7.58)

dove k è l’indice che corre sulle n(n+1)/2 coppie i(k), j(k) di funzioni gaussiane. La matrice di overlap Sekk0 pu` o essere scritta in termini delle Sij per il caso idrogenoide definite dall’Eq.(6.90):

Sekk0 = hBk |Bk0 i = Sii0 Sjj 0 + Sij 0 Sji0 .

(7.59)

e 0 , dell’Hamiltoniano: Gli elementi di matrice, H kk e 0 = hBk |H|B 0 i, H kk k 1 2

H=−

¯h2 ∇21 Zqe2 ¯h2 ∇22 Zqe2 q2 − − − + e 2me r1 2me r2 r12

(7.60)

http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/F90/gauss/helium gauss.f90 http://www.fisica.uniud.it/˜giannozz/Corsi/MQ/Software/C/helium gauss.c

88

si possono scrivere in termini degli elementi di matrice Hij = HijK +HijV ottenuti per il caso idrogenoide con Z = 2, Eq.(6.91) e (6.92):

e 0 = H 0 S 0 + H 0 S 0 + S 0 H 0 + H 0 S 0 + hBk |Vee |B 0 i, H k ij ji jj ii ij ji ii jj kk

(7.61)

e dell’elemento di matrice dell’interazione coulombiana Vee : hBk |Vee |B i = k0

qe2 b 0 (r1 )bj(k0 ) (r2 )d3 r1 d3 r2 (7.62) r12 i(k ) Z q2 + bi(k) (r1 )bj(k) (r2 ) e bj(k0 ) (r1 )bi(k0 ) (r2 )d3 r1 d3 r2 . r12

Z

bi(k) (r1 )bj(k) (r2 )

Ognuno dei due termini in Vee ha la forma Z

I=

2

1 3 3 d r1 d r2 . r12

2

e−αr1 e−βr2

(7.63)

Cerchiamo un cambiamento di variabile che faccia apparire (r1 −r2 )2 all’esponente delle gaussiane: h

αr12 + βr22 = γ (r1 − r2 )2 + (ar1 + br2 )2 



αβ  2 (r1 − r2 ) +  α+β

=

r

i

α r1 + β

(7.64) s

2  β   r2  .

α

(7.65)

Facciamo ora un cambio di variabile da (r1 , r2 ) a (r, s), dove r

r = r1 − r2 ,

s=

α r1 + β

s

β r2 ; α

(7.66)

l’integrale diventa Z

I=

αβ 2 αβ 2 1 r s − α+β − α+β

e

r

e

∂(r1 , r2 ) 3 3 ∂(r, s) d rd s,

(7.67)

dove lo Jacobiano si calcola facilmente come determinante della matrice di trasformazione, Eq.(7.66): ∂(r1 , r2 ) ∂(r, s) =

√

αβ α+β

!3

.

(7.68)

Il calcolo degli integrali è banale e fornisce il risultato cercato: hBk |Vee |Bk0 i =

qe2 π 5/2 qe2 π 5/2 + , αβ(α + β)1/2 α0 β 0 (α0 + β 0 )1/2

(7.69)

dove α = αi(k) +αi(k0 ) , β = αj(k) +αj(k0 ) , α0 = αi(k) +αj(k0 ) , β 0 = αj(k) +αi(k0 ) . (7.70)

89

In modo analogo si calcolano gli elementi di matrice fra funzioni prodotto simmetrizzate, Bk , di gaussiane di tipo P (come quelle definite in Eq.6.88). La combinazione di gaussiane P che ha L = 0 ha la forma: 1 Bk (r1 , r2 ) = √ (r1 · r2 ) bi(k) (r1 )bj(k) (r2 ) + bi(k) (r2 )bj(k) (r1 ) 2

(7.71)

E’ facile verificare che il prodotto di una gaussiana S e una P non contribuisce, in quanto dispari, allo stato fondamentale. Nel caso in cui si usino solo gaussiane di tipo S, il codice scrive su file (”gs-wfc.out”) la funzione: P (r1 , r2 ) = (4πr1 r2 )2 |ψ(r1 , r2 )|2 ,

(7.72)

dove P (r1 , r2 )dr1 dr2 è la probabilità congiunta di trovare un elettrone fra r1 e r1 + dr1 e un elettrone fra r2 e r2 + dr2 . La probabilità di trovare un elettrone fra r e r + dr è data da p(r)dr, con p(r) = 4πr2

Z

|ψ(r, r2 )|2 4πr22 dr2 =

Z

P (r, r2 )dr2 .

(7.73)

E’ facile vedere che per una funziona d’onda prodotto di due funzioni identiche, come la 7.47 o in generale la 7.35, P (r1 , r2 ) = p(r1 )p(r2 ).

7.10.1

Laboratorio

• osservare gli effetti legati al numero di funzioni di base, alla scelta dei coefficienti delle gaussiane λ, all’inclusione di gaussiane P • confrontare l’energia con quella ottenuta con altri metodi: teoria perturbativa con funzioni d’onda idrogenoidi (sezione 7.8), teoria variazionale con Z effettivo (sezione 7.9), risultato esatto (-5.8074 Ry). • fate un plot della P (r1 , r2 ) e della differenza P (r1 , r2 )−p(r1 )p(r2 ), usando per esempio gnuplot e i seguenti comandi: set view 0, 90 unset surface set contour set cntrparam levels auto 10 splot [0:4][0:4] "gs-wfc.out" u 1:2:3 w l Notate come la probabilità P (r1 , r2 ) (la colonna 3 in ”splot”) non sia esattamente uguale al prodotto di p(r1 ) e p(r2 ) (colonna 4; la colonna 5 è la differenza fra i due).

90

Chapter 8

Metodo di Hartree-Fock 8.1

Il metodo di Hartree

L’idea del metodo di Hartree è quella di provare ad approssimare “nel modo migliore possibile” la funzione d’onda soluzione dell’equazione di Schrödinger a N elettroni (problema la cui soluzione esatta è in pratica impossibile da ottenere) con un prodotto di funzioni d’onda a un elettrone—che chiameremo orbitali atomici. Come abbiamo visto, il modo migliore possibile consiste nell’applicare il principio variazionale, minimizzando l’energia media dello stato hψ|H|ψi.

8.1.1

Definizioni

L’hamiltoniano di un atomo con nucleo di carica Z ed N elettroni è H=−

X ¯ h2 i

2me

∇2i −

X Zq 2 e i

ri

+

X q2 e

hiji

rij

(8.1)

dove le parentesi angolari in hiji significano che la somma va limitata a tutte le coppie, ossia implica che i 6= j e che ciascuna coppia viene considerata una volta sola. Ad esempio una scelta valida potrebbe essere X

=

N −1 X

N X

(8.2)

i=1 j=i+1

hiji

Per futura economia di notazione, e anche per maggiore generalità delle formule che ricaveremo, definiamo gli operatori a uno e a due elettroni come fi ≡ − gij

≡

¯h2 2 Zqe2 ∇ − 2me i ri

qe2 rij

(8.3) (8.4)

Con questa notazione l’hamiltoniano si scrive H=

X

fi +

i

X hiji

91

gij

(8.5)

Per usare questa notazione in modo generale è necessario assumere che g agisca simmetricamente sui due elettroni (cosa senz’altro vera per l’interazione coulombiana).

8.1.2

Equazioni di Hartree

Supponiamo allora che la funzione d’onda totale sia esprimibile in forma di prodotto di funzioni d’onda a un elettrone, che supporremo ortonormali (e quindi diverse fra di loro): ψ(1, 2, . . . , N ) = φ1 (1)φ2 (2) . . . φN (N )

(8.6)

Z

φi (1)φj (1) dv1 = δij .

(8.7)

In φi (i), l’argomento i significa: le variabili associate all’elettrone i, ossia in pratica la sua posizione e la proiezione del suo momento angolare di spin. dvi indica un’integrazione rispetto a tutte queste variabili (e quindi anche una somma sugli spin). Invece l’indice i significa: l’insieme dei numeri quantici con cui cataloghiamo una particolare autofunzione, ossia in pratica (nel caso del campo centrale, che come vedremo è quello che ci interessa) il numero quantico principale n, il momento angolare orbitale ` e la sua proiezione m. Procediamo col calcolare il valor medio dell’energia: hψ|H|ψi =

  X X φ∗1 (1) . . . φ∗N (N )  fi + gij  φ1 (1) . . . φN (N ) dv1 . . . dvN

Z

i

=

XZ

φ∗i (i)fi φi (i) dvi +

i

=

XZ

hiji

XZ

φ∗i (i)φ∗j (j)gij φi (i)φj (j) dvi dvj

hiji

φ∗i (1)f1 φi (1) dv1 +

i

XZ

φ∗i (1)φ∗j (2)g12 φi (1)φj (2) dv1 dv2

hiji

(8.8) Nel primo passaggio si è fatto uso dell’ortonormalità (8.7); nel secondo si è semplicemente effettuata una sostituzione delle variabili di integrazione, che non giocano alcun ruolo, standardizzandole a 1 e 2 per comodità. Applichiamo ora il principio variazionale nella formulazione (6.45), richiedendo che nella variazione restino stazionari anche tutti gli integrali Z

Ik =

φ∗k (1)φk (1) dv1

(8.9)

affinchè la normalizzazione di ogni funzione orbitale venga preservata. Vogliamo quindi imporre ! δ hψ|H|ψi −

X

k Ik

=0

(8.10)

k

dove gli k sono i moltiplicatori di Lagrange da determinare. Supponiamo che venga fatta variare soltanto la funzione orbitale con indice k. Si ha allora Z

δIk =

δφ∗k (1)φk (1) dv1 + c.c.

92

(8.11)

(mentre le variazioni di tutti gli altri integrali di normalizzazione saranno ovviamente nulle) ed anche, seguendo le argomentazioni in 6.3.1 sull’hermiticità di H, Z

δhψ|H|ψi = +

δφ∗k (1)f1 φk (1) dv1 + c.c.

XZ

δφ∗k (1)φ∗j (2)g12 φk (1)φj (2) dv1 dv2 + c.c.

(8.12)

j6=k

Questo si ottiene notando che gli unici termini della (8.8) interessati sono quelli con i = k oppure j = k, e che ogni coppia è contata una e una sola volta. Ad esempio per 4 elettroni le coppie sono 12, 13, 14, 23, 24, 34; se scelgo k = 3 P viene un contributo solo da 13, 23 e 34, che corrisponde a una j6=k (l’ordine in cui appaiono gli indici della coppia non ha importanza, essendo g un operatore simmetrico). Pertanto, il principio variazionale assume la forma  Z



δφ∗k (1) f1 φk (1) +

XZ

φ∗j (2)g12 φk (1)φj (2) dv2 − k φk (1) dv1 + c.c. = 0

j6=k

Affinchè il principio variazionale sia soddisfatto, δφk deve essere considerata una variazione arbitraria, e quindi il termine tra le parentesi quadrate deve essere nullo (assieme al suo complesso coniugato!). Deve quindi essere verificata XZ

f1 φk (1) +

φ∗j (2)g12 φk (1)φj (2) dv2 = k φk (1)

(8.13)

j6=k

` utile scriverle metQueste (per k = 1, . . . , N ) sono le equazioni di Hartree. E tendo gli operatori in forma esplicita: 2

−

Zqe2





¯h ∇2 φk (1) − φk (1) +  2me 1 r1 j6=k

XZ

φ∗j (2)

2 φj (2) dv2  φk (1) = k φk (1) r12

(8.14) Osserviamo che ciascuna di esse è simile ad una equazione di Schrödinger, in cui al potenziale coulombiano si aggiunge il “potenziale di Hartree” Z

VH (r1 ) =

ρk (2)

qe2 dv2 r12

(8.15)

dove si è posto ρk (2) =

X

φ∗j (2)φj (2)

(8.16)

j6=k

a dovuta a tutti gli elettroni diversi da quello per cui stiamo ρj è la densit` scrivendo l’equazione.

93

8.1.3

Significato del potenziale di Hartree

La (8.15) rappresenta il potenziale elettrostatico nel punto r1 generato da una distribuzione spaziale di carica ρk . E’ chiaro dunque il significato dell’approssimazione di Hartree. Assumendo che la ψ sia fattorizzabile in un prodotto, abbiamo in pratica assunto che gli elettroni siano fra loro indipendenti da un punto di vista formale. Gli elettroni naturalmente non sono affatto indipendenti e interagiscono fra loro in modo intenso. L’approssimazione fatta non è però cos`ı cattiva se si fa “rientrare” la repulsione coulombiana fra elettroni sotto forma di un campo medio VH fra l’elettrone e il nucleo. VH contiene dunque l’effetto combinato di repulsione dell’elettrone che stiamo considerando da parte di tutti gli altri elettroni. Questo effetto si somma all’attrazione coulombiana esercitata dal nucleo, e la scherma parzialmente. Dunque è circa come se gli elettroni fossero indipendenti, ma interagissero tra loro attraverso il potenziale −Zqe2 /r + VH (r) anzichè il solo −Zqe2 /r.

8.1.4

Campo autoconsistente

VH (r) non è un “vero” potenziale, in quanto la sua definizione dipende dalla distribuzione di densit` a degli elettroni, che a sua volta dipende dalle funzioni orbitali soluzioni della nostra equazione. Il potenziale dunque non è noto a priori, ma è funzione della soluzione; questo tipo di equazione è noto come equazione integro-differenziale. L’equazione pu` o essere risolta in modo iterativo, dopo aver assunto una condizione iniziale per le funzioni orbitali. La procedura è la seguente: 1. calcolare la densit` a di carica (somma dei moduli quadri delle funzioni d’onda) 2. calcolare il potenziale di Hartree generato da questa densità (con metodi di elettrostatica) 3. risolvere l’equazione di Schrödinger ottenendo le funzioni d’onda. La risoluzione dell’equazione può essere effettuata utilizzando i metodi presentati nel capitolo 5.1. La densit` a elettronica è costruita popolando le funzioni d’onda in ordine di energia crescente (soddisfacendo al principio di Pauli!) fino a che tutti gli elettroni sono stati sistemati. In generale le funzioni d’onda ottenute alla fine della procedura sono diverse da quelle iniziali. La procedura può però essere iterata (sia direttamente, usando come funzioni iniziali quelle finali dell’iterazione precedente, che con metodi pi` u raffinati) e porta ad una convergenza del risultato. E’ dunque possibile ripetere la procedura descritta finchè tutte le quantità sono sostanzialmente identiche a quelle dell’iterazione precedente. Il campo VH risultante è allora consistente con le funzioni d’onda, e per questo motivo questo viene chiamato metodo del campo autoconsistente. In calcoli atomici spesso esiste un’ulteriore semplificazione: VH è un campo centrale, ossia dipende dalla sola distanza r1 dell’elettrone dal nucleo. Quando ciò non è vero, si pu` o comunque imporre come approssimazione: basta prendere

94

la media sferica nella definizione di ρk . Se le funzioni d’onda sono soluzione di un problema in un campo centrale, sappiamo a priori che saranno fattorizzate secondo la (5.45). La parte angolare è costituita dalle armoniche sferiche, identificate dai numeri quantici ` e m, mentre la parte radiale è caratterizzata dai numeri quantici n e `. Ovviamente non c’è pi` u la degenerazione dell’energia per diversi `.

8.1.5

Autovalori ed Energia di Hartree

Moltiplichiamo l’equazione di Hartree, Eq(8.13), per φ∗k (1), integriamo e sommiamo: si ottiene X k

k =

XZ

φ∗k (1)f1 φk (1)dv1 +

k

XXZ

φ∗k (1)φk (1)g12 φ∗j (2)φj (2)dv1 dv2 .

k j6=k

(8.17) Confrontiamo questa espressione con l’energia del sistema di molti elettroni, Eq.(8.8). L’energia di repulsione coulombiana è contata due volte, in quanto ogni coppia < jk > appare due volte nella somma. L’energia è quindi data dalla somma degli autovalori dell’equazione di Hartree, meno l’energia di repulsione coulombiana: E=

X k

8.2

k −

X Z

φ∗k (1)φk (1)g12 φ∗j (2)φj (2)dv1 dv2 .

(8.18)

Il metodo di Hartree-Fock

Il metodo di Hartree, che si basa sull’assunzione (8.6), costruisce una funzione d’onda totale che non ha la proprietà di essere antisimmetrica per scambio di una coppia. In base a quanto discusso all’inizio di questo capitolo, è evidente come sia pi` u desiderabile lavorare con una forma antisimmetrica, ossia con un determinante di Slater: φ1 (1) . . . φ1 (N ) 1 . . . ψ(1, . . . , N ) = √ . . . N ! φN (1) . . . φN (N )

(8.19)

La variante del metodo che utilizza funzioni di questo genere è quella comunemente usata, ed è nota come metodo di Hartree-Fock. Si pu` o ripercorrere per la (8.19) tutto lo schema seguito per arrivare alle equazioni di Hartree (8.13). Le complicazioni sono puramente algebriche, legate ` di molto aiuto la proprietà, valida per qualsiasi alla funzione determinantale. E operatore F e funzioni determinantali ψ e ψ 0 : Z φ∗1 (1) . φ∗1 (N ) 1 . . . hψ|F |ψ 0 i = F N ! ∗ φN (1) . φ∗N (N ) φ0 (1) Z 1 ∗ ∗ . = φ1 (1) . . . φN (N )F 0 φN (1)

95

φ0 (1) . φ0 (N ) 1 1 . . . dv1 . . . dvN 0 φN (1) . φ0N (N ) . φ01 (N ) . . (8.20) dv1 . . . dvN 0 . φN (N )

(ossia, sviluppando il primo determinante si ottengono N ! termini che, una volta integrati, sono identici tra loro). Da questa proprietà si ottengono subito i prodotti scalari che ci interessano relativi agli operatori ad uno e due elettroni: hψ|

X

fi |ψi =

i

XZ

φ∗i (1)f1 φi (1) dv1

(8.21)

i

(come nel caso Hartree), e hψ|

X

gij |ψi =

hiji

XZ

φ∗i (1)φ∗j (2)g12 [φi (1)φj (2) − φj (1)φi (2)] dv1 dv2

(8.22)

hiji

le cui dimostrazioni sono semplici anche se noiose. Negli integrali si intende sempre inclusa anche la somma sugli spin. Supponiamo ora che in quest’ultimo termine l’operatore g12 dipenda solo dalle coordinate, come nel caso dell’interazione coulombiana. In tal caso il secondo termine Z

φ∗i (1)φ∗j (2)g12 φj (1)φi (2) dv1 dv2

(8.23)

deve essere nullo nel caso in cui gli spin degli stati i e j siano diversi. Infatti, pensiamo di fattorizzare le funzioni d’onda in una parte dipendente dalle sole coordinate e una parte di spin e notiamo che g12 non modifica quest’ultima. La parte di spin dell’integrale è allora il prodotto di due prodotti scalari fra stati di spin diverso, che per l’ortonormalità sono nulli. Passando ad uno schema in cui invece le variabili di spin non sono esplicitamente incluse, la (8.22) si pu` o scrivere hψ|

X

gij |ψi =

hiji

XZ

φ∗i (1)φ∗j (2)g12 [φi (1)φj (2) − δ(σi , σj )φj (1)φi (2)] dv1 dv2

hiji

(8.24) dove σi è lo spin dell’elettrone i, e δ(σi , σj ) = 0 se σi 6= σj = 1 se σi = σj Riassumendo: hψ|H|ψi =

XZ

φ∗i (1)f1 φi (1) dv1

(8.25)

i

+

XZ

φ∗i (1)φ∗j (2)g12 [φi (1)φj (2) − δ(σi , σj )φj (1)φi (2)] dv1 dv2

hiji

Si procede quindi all’applicazione del principio variazionale. A rigore è necessario imporre non solo che tutte le φi restino normalizzate a 1, ma anche che tutte le coppie φi , φj con lo stesso spin restino fra loro ortogonali. Quest’ultima condizione non era necessaria per il metodo di Hartree e per questo motivo non l’abbiamo menzionata prima. Questo genera una matrice (triangolare) di moltiplicatori di Lagrange ij . Tuttavia, si può far vedere (i dettagli sono sul libro di J. C. Slater, Teoria quantistica della materia, Zanichelli, 1980). che si può

96

sempre pensare di prendere una soluzione per cui la matrice degli è diagonale mediante una semplice trasformazione. Supporremo che sia stata fatta questa scelta. Omettiamo i dettagli e diamo direttamente le risultanti equazioni di HartreeFock, ottenute al solito pensando di variare una sola funzione φk : f1 φk (1) +

XZ

φ∗j (2)g12 [φk (1)φj (2) − δ(σk , σj )φj (1)φk (2)] dv2 = k φk (1)

j

(8.26) o in forma esplicita −

XZ ¯h2 2 Zqe2 q2 ∇1 φk (1) − φk (1) + φ∗j (2) e [φj (2)φk (1) 2me r1 r12 j

−

(8.27)

δ(σk , σj )φk (2)φj (1)] dv2 = k φk (1)

L’energia del sistema, Eq. 8.25, si può esprimere, in modo analogo a come visto per il caso di Hartree, tramite la somma degli autovalori della 8.27, meno un termine che compensa il doppio conteggio dell’energia di repulsione coulombiana e dell’energia di scambio: E=

X

k −

k

X Z

φ∗k (1)φ∗j (2)g12 [φk (1)φj (2) − δ(σj , σk )φj (1)φk (2)] dv1 dv2 .

(8.28) Osserviamo attentamente le differenze rispetto alle equazioni di Hartree (8.13): 1. la

P

j

comprende anche il caso j = k.

2. per gli elettroni j con lo spin identico a quello di k c’è un termine in pi` u, detto termine di scambio 3. a causa del termine di scambio, il caso j = k dà un contributo non nullo solo se gli spin sono diversi. Cerchiamo di capire cosa c’è dietro. Prima di procedere osserviamo che la (8.27) avr` a normalmente infinite soluzioni, di cui solo le N a energia pi` u bassa verranno occupate da elettroni. Gli stati che restano liberi sono gli stati eccitati. P La j va pensata limitata agli stati occupati.

8.2.1

Potenziale colombiano e di scambio

Riscriviamo l’equazione di Hartree-Fock sotto la forma −

¯h2 2 Zqe2 ∇1 φk (1) − φk (1) + VH (1)φk (1) + (Vˆx φk )(1) = k φk (1), 2me r1

(8.29)

dove abbiamo definito un potenziale coulombiano (o anche ”di Hartree”, ma non è lo stesso che nel metodo di Hartree!) VH ed un potenziale di scambio Vx . Il potenziale coulombiano è lo stesso per tutti gli orbitali: VH (1) =

XZ j

φ∗j (2)

qe2 φj (2)dv2 ≡ r12

97

Z

ρ(2)

qe2 dv2 , r12

(8.30)

dove abbiamo introdotto la densità di carica: ρ(2) =

X

φ∗j (2)φj (2).

(8.31)

j

Si pu` o verificare che ρ è uguale alla probabilità di trovare un elettrone: Z

ρ(1) = N

|Ψ(1, 2, .., N )|2 dv2 ...dvN .

(8.32)

Il termine di scambio: (Vˆx φk )(1) = −

X

Z

δ(σk , σj )

φ∗j (2)

j

qe2 φk (2)φj (1) dv2 r12

(8.33)

non ha la forma semplice del potenziale coulombiano: VH (1)φk (1), dove VH (1) contiene l’integrazione sulla variabile 2. Ha piuttosto una forma del tipo (Vˆx φ)(1) ≡

Z

Vx (1, 2)φk (2)dv2

(8.34)

che caratterizza una interazione non locale.

8.2.2

La densit` a di scambio

Per cercare di comprendere il termine di scambio, definiamo la seguente “densità di scambio”: X φ∗ (1)φ∗j (2)φj (1)φk (2) ρx (2) ≡ δ(σk , σj ) k (8.35) φ∗k (1)φk (1) j Con questa definizione, il termine di scambio (8.33) si può riscrivere (Vˆx φ)(1) = −

"Z

#

q2 ρx (2) e dv2 φk (1) r12

(8.36)

Lo scopo di questa definizione un po’ artificiosa è stato quello di riuscire a scrivere formalmente il termine di scambio come il prodotto di un “potenziale efficace” per φk (1), risultante da interazioni elettrostatiche con una distribuzione spaziale di densit` a di carica, esattamente come abbiamo fatto per il termine di Hartree. La densit` a di scambio ha le seguenti proprietà: 1. vi contribuiscono solo gli elettroni con lo stesso spin di quello che stiamo considerando 2. rappresenta una quantit` a di carica totale pari a 1. Infatti Z

ρx (2)dv2 =

X j

φ∗ (1)φj (1) δ(σk , σj ) ∗k φk (1)φk (1)

Z

φ∗j (2)φk (2) dv2

(8.37)

da cui si vede che i termini con j 6= k danno contributo nullo, mentre quello con j = k d` a un contributo pari a uno. Questa carica è dunque quella dell’elettrone che sto considerando, ma distribuita in qualche modo nello spazio. Da notare che se considero uno stato k eccitato, ossia non occupato, allora è sempre j 6= k e la carica di scambio è nulla, come si conviene ad un elettrone . . . che non c’è.

98

3. esaminiamo cosa succede quando i punti 1 e 2 tendono a coincidere. In questo limite la carica di scambio tende a ρσk (1) =

X

δ(σk , σj )φ∗j (1)φj (1)

(8.38)

j

ossia alla densit` a totale nel punto 1 di tutti gli elettroni con lo stesso spin dell’elettrone k. Vediamo dunque che, nell’equazione di Hartree-Fock per l’elettrone k, il primo termine (quello di Hartree, ma con somma estesa a tutti i j) include le interazioni efficaci con tutti gli elettroni, incluso k. Il termine di scambio toglie l’interazione dell’elettrone che stiamo considerando con se stesso. Consideriamo ora la densit` a di carica totale (Hartree+scambio) con lo stesso spin dell’elettrone k. Questa densità è nulla nel punto in cui si trova k, perchè in quel punto il termine di scambio compensa esattamente il termine di Hartree come indicato dall’ultima proprietà sopra. Quindi è come se il nostro elettrone trascinasse con sè una buca, detta buca di Fermi, che tiene lontani gli altri elettroni con lo stesso spin. Si tratta ovviamente di un effetto dovuto al principio di esclusione. Notiamo infine che il metodo di Hartree il termine escluso dalla somma j = k esclude—in modo pi` u rudimentale—l’interazione dell’elettrone k con se stesso: lo stesso effetto rappresentato dal termine di scambio nelle equazioni di Hartree-Fock.

8.2.3

L’atomo di elio

Come gi` a visto per il metodo di Hartree, anche il metodo di Hartree-Fock utilizzato comunemente per gli atomi adotta l’approssimazione di campo centrale. Ciò permette la fattorizzazione delle Eq.(8.27) in una parte radiale e una parte angolare, e la classificazione delle soluzioni con i ”tradizionali” numeri quantici n, `, m. Consideriamo il caso pi` u semplice dell’atomo a due elettroni e confrontiamo le equazioni di Hartree (8.14) con quelle di Hartree-Fock (8.27). L’equazione di Hartree (8.14) si riduce a (considerando k = 1) Zqe2 ¯h2 2 ∇1 φ1 (1) − φ1 (1) + − 2me r1

"Z

#

φ∗2 (2)

qe2 φ2 (2) dv2 φ1 (1) = 1 φ1 (1) (8.39) r12

L’equazione di Hartree-Fock (8.27) invece si riduce a −

¯h2 2 Zqe2 q2 ∇1 φ1 (1) − φ1 (1) + φ∗1 (2) e [φ1 (2)φ1 (1) − φ1 (2)φ1 (1)] dv2 2me r1 r12 Z 2 q + φ∗2 (2) e [φ2 (2)φ1 (1) − δ(σ1 , σ2 )φ1 (2)φ2 (1)] dv2 = 1 φ1 (1) (8.40) r12 Z

ovvero, poichè l’integrando nel primo integrale è nullo, ¯h2 2 Zqe2 q2 − ∇1 φ1 (1) − φ1 (1) + φ∗2 (2) e [φ2 (2)φ1 (1) 2me r1 r12 −δ(σ1 , σ2 )φ1 (2)φ2 (1)] dv2 = 1 φ1 (1). Z

99

(8.41)

Assumiamo ora di andare alla ricerca dello stato fondamentale. In base alle considerazioni effettuate nella sezione 7.7, ci aspettiamo che nello stato fondamentale i due elettroni abbiano spin opposto (δ(σ1 , σ2 ) = 0) e che entrambi occupino lo stesso orbitale 1s (ossia che φ1 e φ2 siano la stessa funzione φ1 ). Questo fa s`ı che l’equazione di Hartree-Fock (8.41) per lo stato fondamentale sia identica a quella di Hartree (8.39): i metodi di Hartree ed Hartree-Fock sono quindi in questo caso equivalenti. Questo avviene in sostanza perchè i due elettroni hanno spin opposto e quindi non vi sono effetti di scambio in azione. In generale, si parla di Hartree-Fock Ristretto (RHF) per il caso frequente in cui tutti gli orbitali sono presenti a coppie, formate da una stessa funzione di r moltiplicata per funzioni di spin opposto.

8.3

L’energia di correlazione

La soluzione Hartree-Fock non è esatta: lo sarebbe se il sistema sotto esame fosse descritto da una funzione d’onda formata da un solo determinante di Slater. In generale, ci` o non è vero. La differenza di energia fra la soluzione esatta e la soluzione Hartree-Fock va sotto il nome di energia di correlazione.1 L’origine del nome deriva dal fatto che nell’approssimazione di Hartree-Fock manca una parte della ”correlazione elettronica”, ovvero dell’effetto che un elettrone ha sugli altri. Quest’ultimo è presente tramite l’interazione di scambio e l’interazione elettrostatica, ma mancano effetti pi` u sottili che si riflettono in una forma della funzione d’onda esatta che è pi` u generale di quella Hartree-Fock. Abbiamo visto tali effetti in opera nel caso dell’atomo di He risolto con il codice helium gauss (sezione 7.10): la probabilit` a P (r1 , r2 ) di trovare un elettrone a distanza r1 e uno a distanza r2 dal centro non è semplicemente uguale a p(r1 )p(r2 ), perchè gli elettroni cercano di ”evitarsi”. L’energia di correlazione nel caso dell’atomo di He è circa 0.084 Ry: una quantità piccola rispetto alle energie in gioco (∼ 1.5% dell’energia) ma non trascurabile. E’ lecito chiedersi a cosa possa mai servire una simile quantità, il cui calcolo è possibile solo in sistemi molto semplici (come per esempio l’atomo di Elio) o in sistemi modello (per esempio: il gas omogeneo di elettroni) nei quali si riesce a ottenere una soluzione praticamente esatta. Tuttavia avere una stima di quanto si sbaglia e soprattutto in che casi si sbaglia di pi` u e perchè è cosa assai utile, nonchè la base di partenza per migliorare il metodo. Un ovvio modo per migliorare i risultati Hartree-Fock è aggiungere contributi da altri determinanti di Slater alla funzione d’onda di prova. Questa è l’essenza del metodo della ”interazione delle configurazioni”. La sua applicazione pratica richiede una raffinata ”tecnologia” per scegliere fra l’enorme numero di possibili determinanti di Slater il sottoinsieme pi` u significativo. E’ molto usato in chimica quantistica per ottenere risultati di alta precisione, ma si tratta di un approccio computazionalmente molto pesante edi fatto limitato a sistemi (molecole) semplici. Metodi pi` u agili danno una stima abbastanza buona dell’energia di correlazione tramite la teoria delle perturbazione (i cosid1

Feynman la chiamava energia di stupidit` a, perché la sola grandezza fisica che misura è la nostra incapacit` a di trovare la soluzione esatta!

100

detti approcci MP, Møller-Plesset). Un approccio completamente diverso è invece quello della teoria del funzionale densità, benché alla fine produca delle equazioni molto reminiscenti delle equazioni Hartree-Fock.

8.4

Programma: helium hf radial

Il programma helium hf radial.f902 (oppure helium hf radial.c3 ) risolve le equazioni di Hartree-Fock (equivalenti in questo caso a quelle di Hartree come visto sopra) per lo stato fondamentale dell’atomo di elio. helium hf radial è basato su hydrogen ed utilizza lo stesso algoritmo di integrazione basato sul metodo di Numerov. La parte nuova è costituita dall’implementazione del metodo del campo autoconsistente per la ricerca degli orbitali. Il calcolo consiste nella risoluzione dell’equazione di Schrödinger radiale sotto un potenziale effettivo Vscf , somma del potenziale coulombiano del nucleo pi` u il potenziale di Hartree come definito sopra: Zq 2 Vscf (r) = − e + VH (r), r

q2 VH (r1 ) = e 2

Z

|R(r1 )|2 3 d r2 . r12

(8.42) (0)

Si parte da una stima iniziale di VH (r) nella routine init pot ( VH (r) = 0, semplicemente). Con lo stato fondamentale R(r) ottenuto da tale potenziale, si calcola (nella routine rho of r) la densità di carica ρ(r) = |R(r)|2 . La routine v of rho ricalcola il nuovo potenziale di Hartree VHout (r) per semplice integrazione, usando il teorema di Gauss: VeHout (r) = V0 + qe2

Z r rmax

Q(s) ds, s2

Z

Q(s) =

ρ(r)4πr2 dr

(8.43)

r