Post on 25-Feb-2021
Facolta di Scienze Matematiche, Fisiche e Naturali
Laurea Triennale in Fisica
Effetti relativistici negli spettriatomici: un approccio
computazionale
Relatore: Dott. Nicola Manini
Francesco Brivio
Matricola n◦ 715099
A.A. 2008/2009
Codice PACS: 31.15.aj
ai miei genitori
Effetti relativistici negli spettri
atomici: un approccio
computazionale
Francesco Brivio
Dipartimento di Fisica, Universita degli Studi di Milano,
Via Celoria 16, 20133 Milano, Italia
11 Dicembre 2009
Sommario
Calcoliamo numericamente gli autovalori di energia dell’equazione di
Dirac per un singolo elettrone per un qualsiasi potenziale centrale V (r).
Cio ci permette di valutare quantitativamente gli effetti relativistici nel
moto degli elettroni degli atomi. Il codice da noi sviluppato ottiene gli
autovalori ed autofunzioni di Dirac integrando un’equazione differenziale
ordinaria cui sono imposte opportune condizioni al contorno.
Relatore: Dr. Nicola Manini
5
Indice
1 Introduzione 7
2 Teoria 9
2.1 L’equazione di Dirac . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Derivazione dell’equazione radiale . . . . . . . . . . . . . . . . . . 9
2.3 Derivazione dell’equazione da implementare numericamente . . . . 12
3 Implementazione numerica 14
3.1 L’algoritmo Shooting . . . . . . . . . . . . . . . . . . . . . . . . . 15
4 Risultati 17
4.1 Potenziale Coulombiano a nucleo puntiforme . . . . . . . . . . . . 17
4.2 Potenziale Coulombiano a nucleo finito . . . . . . . . . . . . . . . 21
4.3 Effetti relativistici in presenza di screening . . . . . . . . . . . . . 22
5 Discussioni e Conclusioni 28
Bibliografia 29
6
1 Introduzione
Per determinare lo spettro dei livelli elettronici in un atomo, normalmente, si
risolve un’equazione agli autovalori del tipo
H0ψ = E0ψ H0 =p
2
2m+ V (r) (1)
con un opportuno potenziale centrale V (r), che per atomi ad un elettrone e pari
al potenziale Coulombiano −Ze2
r. Risolvendo in coordinate polari [1] si separa la
parte angolare dalla parte radiale. Per la parte radiale si ottiene cosı l’equazione
di Schrodinger[−
~2
2m▽2 +
~2l(l + 1)
2mr2+ V (r)
]Rnl = ERnl (2)
in cui considero la massa del nucleo mn ≫ m (massa dell’elettrone) per cui la
massa ridotta µ ≃ m e la componente radiale dell’operatore Laplaciano ▽2 in
coordinate polari sferiche definito nel seguente modo
▽2 =1
r2
d
dr
(r2 d
dr
)=
2
r
d
dr+
d
dr2.
Risolvendo analiticamente l’equazione si ottengono i valori di energie per-
messe [1]
E0 = εn = −m(Ze2)2
2~2n2. (3)
Questo approccio e valido nel limite non relativistico in cui la velocita dell’elet-
trone e trascurabile rispetto a quella dalla luce c. Di fatto pero gli elettroni
negli atomi, nelle regioni vicino al nucleo, raggiungono velocita non lontane da c.
Risulta quindi necessario considerare le correzioni relativistiche all’Hamiltoniana
classica considerata nel problema di Schrodinger di equazione (1).
Gli effetti relativistici negli spettri atomici sono piccoli, ma crescono ra-
pidamente con il numero atomico Z, come Z4 [2]. Le correzioni allo spettro,
soprattutto per atomi pesanti (grande Z), non sono quindi trascurabili e diventa
importante prenderle in considerazione.
Un possibile modo di procedere per trattare il problema relativistico e di
utilizzare la teoria delle perturbazioni, considerando il termine relativistico del-
l’Hamiltoniana come una perturbazione al termine classico. L’Hamiltoniana da
considerare e dunque
H = H0 + Hrel , (4)
7
dove H0 e l’Hamiltoniana classica (1) e Hrel contiene i termini della correzione
relativistica all’ordine piu basso (il secondo) in vc. Considerando (1) e utilizzando
la teoria delle perturbazioni l’autovalore di energia e
E = E0 + ∆E ∆E = 〈ψ|Hrel|ψ〉 (5)
dove E0 e l’autovalore dell’Hamiltoniana H0 (1) corrispondente all’autostato |ψ〉,
e ∆E e la correzione relativistica agli autovalori di energia classici.
Un modo alternativo e preferibile di procedere, che seguiremo in questo la-
voro di tesi, e di risolvere l’equazione di Dirac. Questa equazione, formulata da
Paul Dirac [3], rappresenta il moto dell’elettrone in maniera interamente relativis-
tica, e conduce all’equazione di Schrodinger nel limite di c→ ∞. Ci proponiamo
un approccio numerico all’equazione di Dirac, in modo da poter calcolare gli
spettri in presenza di un potenziale radiale qualunque.
Adottiamo il sistema di unita di misura atomico. In questo sistema tutte le
grandezze sono derivate ponendo
~ = 1 m = 1 e2 =q2e
4πε0= 1
dove m e la massa a riposo dell’elettrone, qe la carica, ε0 la permittivita del
vuoto e ~ la costante di Planck. In questo sistema di unita, le lunghezze vengono
misurate in raggi di Bohr a0, l’energia in unita di Hartree EHa e le velocita in
unita della velocita di Bohr vB:
a0 =~
2
me2= 0.52917720859 10−10 m
EHa =e2
a0
=me4
~2= 27.21138386 eV
vB =
√EHa
m=e2
~= 2.187691254 106 m/s.
Si nota che vB e 2 ordini di grandezza piu piccola della velocita della luce. Intro-
ducendo la costante di struttura fine α, che misura l’importanza della correzione
relativistica al moto non relativistico degli elettroni,
α =vB
c= 0.007297353 =
1
137.035999679,
otteniamo c nel sistema di unita atomiche:
c =1
αvB = 137.035999679 .
8
2 Teoria
2.1 L’equazione di Dirac
Consideriamo l’equazione di Dirac stazionaria per un singolo elettrone [4]
Hψ = Wψ (6)
con
H = c~α · ~p+ βmc2 + V (r)I . (7)
In queste equazioni ψ e un vettore a quattro componenti:
ψ =
ψ1
ψ2
ψ3
ψ4
=
(ψA
ψB
)
, ψA =
(ψ1
ψ2
)
, ψB =
(ψ3
ψ4
)
,
~α e β sono matrici 4×4:
~α =
(0 ~σ
~σ 0
)
, β =
(I 0
0 −I
)
,
il simbolo I indica l’identita 2×2 o 4×4 a seconda del contesto, il vettore ~σ
rappresenta le matrici di Pauli ~σ = (σx, σy, σz) e ~p rappresenta l’operatore impulso
(px, py, pz) = −i~▽ = −i~(
∂∂x, ∂
∂y, ∂
∂z
). Ora l’equazione (6) si scrive
c~α · ~p
(ψA
ψB
)=(W − βmc2 − V (r)I
)(ψA
ψB
). (8)
Sostituendo i vari termini, (8) puo essere riscritta formalmente in forma matriciale
c
((~σ · ~p) ψB
(~σ · ~p) ψA
)=
(W − V −mc2 0
0 W − V +mc2
)(ψA
ψB
). (9)
2.2 Derivazione dell’equazione radiale
Considerando un potenziale centrale V (r) e possibile separare la parte angolare
dalla parte radiale dell’equazione. Innanzitutto possiamo riscrivere [4]
~σ · ~p =1
r
~σ · ~x
r
(−i~r
∂
∂r+ i~σ · ~L
), (10)
9
e consideriamo
ψA = gχj3k
ψB = ifχj3−k
(11)
in cui f e g sono funzioni che dipendono solo da r, mentre χj3−k e χj3
k sono spinori
a due componenti che dipendono da (ϑ, ϕ). La scelta di considerare (11) e fatta
in modo tale che ψA e ψB abbiano lo stesso j e j3 e differiscano solamente per l
(che puo assumere i valori l = j ± 12
corrispondenti [4] a ±κ). Sostituendo (11) e
(10) in (9) otteniamo
c1
r
~σ · ~x
r
(−i~r
∂
∂r+ i~σ · ~L
)(ifχj3
−k
gχj3k
)
=
=
(W − V −mc2 0
0 W − V +mc2
)(gχj3
k
ifχj3−k
)
. (12)
Definiamo ora un nuovo operatore
K = β(~σ · ~L+ ~) .
Si puo mostrare [4] che sia ψA sia ψB sono autostati dell’operatore K:
KψA,B = −~κψA,B (13)
con
κ =
{−l − 1 = −(j + 1
2) se j = l + 1
2
l = +(j + 12) se j = l − 1
2
. (14)
In questo modo possiamo determinare il comportamento dell’operatore ~σ · ~L sugli
stati ψA,B
~σ · ~L
(ψA
ψB
)= (βK − ~)
(ψA
ψB
)=
(K − ~ 0
0 −K − ~
)(ψA
ψB
)=
=
(~(−κ− 1)ψA
~(κ− 1)ψB
)(15)
e possiamo scrivere
c1
r
~σ · ~x
r
((−i~r ∂
∂r+ i(κ− 1)~)ifχj3
−k
(−i~r ∂∂r
+ i(−κ− 1)~)gχj3k
)
=
=
(W − V −mc2 0
0 W − V +mc2
)(gχj3
k
ifχj3−k
)
. (16)
10
Ora l’operatore ~σ·~xr
agisce solamente sulla parte di momento angolare dello
stato. In dettaglio~σ · ~x
rχj3
k = −χj3−k , (17)
~σ · ~x
rχj3−k = −χj3
k (18)
cosı
c1
r
((i~r ∂
∂r− i(κ− 1)~)ifχj3
k
(i~r ∂∂r
+ i(κ + 1)~)gχj3−k
)
=
=
(W − V −mc2 0
0 W − V +mc2
)(gχj3
k
ifχj3−k
), (19)
che possiamo riscrivere
~c1
r
((−r ∂
∂r+ (κ− 1))f 0
0 (r ∂∂r
+ (κ+ 1))g
)(χj3
k
iχj3−k
)
=
=
((W − V −mc2)g 0
0 (W − V +mc2)f
)(χj3
k
iχj3−k
)
(20)
e cancellando da entrambe le parti i termini uguali
~c
(−∂f
∂r+ κ−1
rf
∂g∂r
+ κ+1rg
)=
((W − V −mc2)g
(W − V +mc2)f
). (21)
Questa equazione e l’equazione di Dirac radiale scritta in forma matriciale. E
possibile riscrivere (21) come un sistema di 2 equazioni differenziali di primo
ordine non omogenee e accoppiate [4]
g′κ = −κ + 1
rgκ +
(W − V +mc2)
c~fκ (22)
f ′κ =
κ− 1
rfκ −
(W − V −mc2)
c~gκ (23)
in cui fκ e gκ sono le due componenti radiali delle autofunzioni ψA,B, come mostra-
to in (11), e κ rappresenta il momento angolare relativistico che e legato a l e j
dalle seguenti relazioni [5]:
κ(κ + 1) = l(l + 1) (24)
j = |κ| −1
2= l +
1
2δ (25)
δ = −sgn(κ) = −|κ|
κ(26)
l = j −1
2δ = |κ| −
1
2(1 + δ) . (27)
11
Dalla relazione (25) notiamo come κ = 0 sia un valore fisicamente non permesso
in quanto produrrebbe j < 0 che non ha significato fisico.
2.3 Derivazione dell’equazione da implementare numeri-
camente
Definiamo ora E = W − mc2 in modo che E non contenga l’energia della massa
a riposo dell’elettrone, e introduciamo la massa relativistica
M = m+1
2c2(E − V ) (28)
che e di fatto una funzione della posizione e dell’autovalore incognito E. Con
queste sostituzioni le equazioni (22) e (23) diventano
g′κ = −κ + 1
rgκ + 2Mcfκ (29)
f ′κ =
κ− 1
rfκ +
1
c(V −E)gκ . (30)
Introduciamo ora una nuova funzione φκ definita nel modo seguente
φκ =1
2Mg′κ (31)
e dalla (29) otteniamo [4]
fκ =g′κ
2Mc+κ+ 1
r
gκ
2Mc=φκ
c+κ+ 1
c
gκ
2Mr. (32)
Considerando la derivata della funzione (28) M(r)
M ′ = −V ′
2c2(33)
e derivando l’equazione (32), otteniamo
f ′κ =
φ′κ
c+κ+ 1
c
(g′κ
2Mr−
1
r2
gκ
2M−gκ
2r
−V ′
2M2c2
)=
=φ′
κ
c+κ+ 1
crφκ +
κ+ 1
2M2crgκV ′
2c2−
κ+ 1
2Mcr2gκ . (34)
Confrontando (34) con (30) e (32) moltiplicando tutto per c ottengo
φ′κ +
κ+ 1
rφκ +
κ+ 1
2M2rgκV ′
2c2−κ+ 1
2Mr2gκ =
κ− 1
r
(φκ +
κ + 1
2Mrgκ
)+ (V − E)gκ ,
(35)
12
da cui
φ′κ(r) = −
2
rφκ(r) +
[(κ+ 1) + κ2 − 1
2M(r)r2−
κ+ 1
4M2(r)c2rV ′(r) + V (r) −E
]gκ . (36)
Ricordando la definizione di φκ (31) e facile ottenere
φ′κ =
g′′κ2M
−g′κM
′
2M2(37)
e sostituendo (37) nell’equazione (36) ottengo
g′′κ2M
−g′κM
′
2M2= −
2
rφκ(r) +
[(κ+ 1) + κ2 − 1
2M(r)r2−
κ + 1
4M2(r)c2rV ′(r) + V (r) −E
]gκ
(38)
da cui
g′′κ = −
(2
r−M ′
M
)g′κ +
[κ(κ+ 1)
2Mr2−
κ+ 1
4M2c2rV ′ + V − E
]2Mgκ . (39)
Sostituendo l’espressione esplicita (33) di M ′ ricavo un’equazione per la
funzione incognita gk
g′′κ = −
(2
r+
V ′
2Mc2
)g′κ +
[(V − E) +
κ(κ+ 1)
2Mr2−
κ+ 1
4M2c2rV ′
]2Mgκ . (40)
Ora se suppongo di studiare il caso non relativistico, cioe considero il limite
per c→ +∞ ottengo l’equazione
g′′k +2
rg′k −
[(V −E) +
κ(κ+ 1)
2mr2
]2mgk = 0 (41)
che e precisamente l’equazione radiale di Schrodinger (2) a patto di porre κ(κ +
1) = l(l + 1), come previsto da Eq. (24). l rappresenta il momento angolare
orbitale.
Nel caso non relativistico allora, gκ → Rnl, mentre fκ → 0. Possiamo quin-
di identificare gκ con la componente “grande” delle autofunzioni dell’equazione
di Dirac, mentre fκ con quella “piccola”. Nella trattazione finora svolta siamo
passati dal dover risolvere un sistema di due equazioni differenziali accoppiate
di primo ordine (22) e (23) per fκ e gκ, ad un’equazione differenziale di secondo
ordine (40) per gκ. Se si vuole trovare l’espressione di fκ, una volta trovata la
soluzione gκ(r) di (40) la si sostituisce in Eq. (30) e si risolve l’equazione dif-
ferenziale. Una volta calcolate le funzioni gκ e fκ possiamo ottenere la densita di
probabilita radiale ρ(r) = r2(|fκ|2 + |gκ|
2).
13
3 Implementazione numerica
Ci proponiamo di risolvere numericamente l’equazione differenziale (40) fissate
opportune condizioni al contorno. Per studiare le condizioni al contorno considero
l’andamento dell’autofunzione gκ e della sua derivata [2].
Nell’origine abbiamo
gκ(r) ∼ rlwκ(r) (42)
g′κ ∼ lrl−1wκ + rlw′κ . (43)
I valori di gκ e della sua derivata prima, nell’origine r1 = 0 dipendono quindi dal
momento angolare l
gκ(0) =
{0 se l 6= 0
wκ(0) se l = 0(44)
e
g′κ(0) =
w′κ(0) se l = 0
wκ(0) se l = 1
0 se l 6= 0, 1
(45)
mentre all’infinito gκ → 0. Chiamando
A =(2
r+
V ′
2Mc2
)e B = 2M
[(V −E) +
κ(κ + 1)
2Mr2−
κ+ 1
4M2c2rV ′]
,
l’equazione (40) si scrive
g′′κ = −Ag′κ +Bgκ . (46)
Risolviamo l’equazione (46) mappandola su un problema di prim’ordine nel doppio
delle variabili. Chiamando y(0) = gκ e y(1) = g′κ, (46) puo essere scritta nella
forma equivalente:
y′(0) = y(1)
y′(1) = −Ay(1) +By(0)
y′(2) = 0
(47)
dove ho un sistema di 3 equazioni perche fissate le condizioni al contorno (46) e
risolvibile solo per certi valori di energia incogniti. Scrivo E come terza compo-
nente della funzione vettoriale y (y(2) = E) per cui la terza equazione del sistema
ci dice che l’energia E e costante. In effetti il sistema ha 3 incognite: le funzioni
y(0) e y(1) e l’autovalore di energia E ed e risolvibile univocamente introducendo
una terza equazione.
Normalmente le condizioni al contorno andrebbero applicate in r1 = 0 e
r2 = ∞. Dovendo pero risolvere il problema in uno spazio finito, e volendo evitare
14
eventuali singolarita nell’origine, di fatto integriamo l’equazione tra r1 ≃ 10−12
a0 e r2 ≃ 50 a0.
3.1 L’algoritmo Shooting
Un generico problema con condizioni al contorno si presenta in questa forma:
cercare la soluzione di un sistema di N equazioni differenziali del primo ordine
che soddisfano n1 condizioni nel punto iniziale x1 e n2 nel punto finale x2. Le
equazioni sono:
dy(i)(x)
dx= g(i)(x, y(0), y(1), ..., y(N−1)) i = 0, 1, ..., N − 1 . (48)
Nel punto iniziale x1 la soluzione deve soddisfare le condizioni:
B1j(x1, y(0), y(1), ..., y(N−1)) = 0 j = 0, ..., n1 − 1 (49)
mentre nel punto finale x2 :
B2k(x2, y(0), y(1), ..., y(N−1)) = 0 k = 0, ..., n2 − 1 . (50)
Per risolvere (48) con un algoritmo di integrazione esplicita 1-dimensionale, ho
bisogno di fissare nel punto iniziale x1 N valori iniziali per le y(i), le quali sono
pero soggette solo a n1 condizioni (49). Ci sono quindi n2 = N −n1 valori iniziali
che possono essere assegnati “liberamente”. Possiamo immaginare [6] che questi
parametri liberi siano le componenti di un vettore V appartenente ad uno spazio
vettoriale di dimensione n2. Assegnando “arbitrariamente” le componenti del
vettore V e applicando le n1 condizioni al contorno fissate, possiamo generare
una particolare y(x1). In altre parole (49) la scrivo come
y(i)(x1) = y(i)(x1, V0, V1, ..., Vn2−1) . (51)
Per un fissato vettore V, siamo di fronte ad un problema di Cauchy con tutte
le N condizioni iniziali specificate che, integrato numericamente ad esempio con
un algoritmo del tipo Runge-Kutta [7], ci permette di ottenere cosı la soluzione
y(x2). Definiamo, nel punto x2, un vettore discrepanza F di dimensione n2 che
indica quanto la soluzione trovata y(x2) differisce rispetto alle n2 condizioni (50)
al punto finale x2. La scelta corretta delle componenti di V, che soddisfano le
condizioni al contorno in x2, sara quella che annulla il vettore F. A questo scopo
utilizziamo il metodo di Newton modificato per renderlo globalmente convergente.
Nel nostro specifico problema (N = 3) fornisce n1 = 2 condizioni nell’origine
x1 = 0 (equazioni (44) e (45)) e n2 = 1 condizioni all’infinito (y(0)(x2) = 0), che
15
Figura 1: Schema del funzionamento dell’algoritmo shooting: fissate le
condizioni al contorno in x1 e partendo da un valore arbitrario per il vettore
V, integriamo numericamente le equazioni (48) dal punto iniziale x1 al
punto finale x2. In generale la soluzione y(x2) non soddisfa le condizioni
al contorno nel punto finale x2, ma con successivi aggiustamenti del valore
del vettore V e successive integrazioni, e possibile imporre la condizione
in x2 (cioe annullare il vettore discrepanza F), determinando cosı il valore
corretto di V.
implementiamo come un punto x2 “molto lontano” da 0 se confrontato con le
dimensioni caratteristiche del problema considerato (il raggio di Bohr). Questi
tre vincoli consentono di determinare l’autovalore E.
A questo punto per applicare l’algoritmo shooting bisogna generare il vettore
V: nel nostro caso esso ha dimensione 1 (n2 = N − n1 = 1) e contiene un’ipotesi
per l’autovalore di energia da trovare (cioe il valore iniziale della funzione costante
y(2)(x)). Sapendo che gli autovalori relativistici non sono lontani da quelli non
relativistici, come primo guess di autovalore di E inserisco un valore dato dalla
formula (3). Ci aspettiamo che dopo alcune iterazioni di Newton, l’algoritmo
restituisca il valore esatto dell’autovalore. In seguito, variando la condizione
iniziale per E, otteniamo l’intero spettro discreto dell’energia.
16
Gli algoritmi utilizzati per risolvere le equazioni differenziali e per risolvere
il sistema lineare sono quelli proposti in Ref. [6]. In particolare abbiamo usato il
metodo di Runge-Kutta di ordine cinque con controllo automatico del passo di
integrazione e il metodo di Newton per la risoluzione dell’equazione (50) associata
al problema agli autovalori con condizioni al contorno.
4 Risultati
4.1 Potenziale Coulombiano a nucleo puntiforme
Effettuiamo ora una prima verifica dell’algoritmo. Consideriamo il caso piu sem-
plice di potenziale V (r): quello Coulombiano prodotto da una carica puntiforme
V (r) = −Ze2
r(52)
e poniamo inizialmente Z = 1: l’atomo di idrogeno.
Dalla Tabella 1 si puo notare come i valori delle energie per n da 1 a 5
si discostano da quello non relativistico (3) per una frazione dell’ordine 10−4 ÷
10−5 dell’autovalore. Dalla tabella si puo notare come gli autovalori di energia
risultino degeneri per l = j ± 12
a dato n e j. Questa e la ben nota degenerazione
“accidentale” caratteristica del potenziale Coulombiano. Otteniamo inoltre gli
splitting attesi di spin orbita, ad esempio tra gli stati 2P3/2 e 2P1/2. Con i valori
di energia riportati in Tabella 1, si puo verificare che lo splitting tra i due stati
considerati risulta E|2,1,3/2 〉−E|2,1,1/2 〉 = 0.000045284 eV = 45.284 µeV in accordo
da quanto ottenuto con la formula relativistica all’ordine α2 [2]:
E = E0 + ∆E = −mZ2e4
2~2n2
[1 +
Z2
n
(2
2j + 1−
3
4n
)α2
]. (53)
Figura 2 mostra lo spettro per n = 5 fissato. Sono visibili la degenerazione
rispetto l e lo splitting rispetto a j in accordo con l’equazione (53).
Confrontiamo ora i risultati ottenuti con il nostro metodo con quelli reperibili
in letteratura scientifica [8]. Gli autovalori disponibili sono ottenuti per l = 0,
κ = −1, j = 1/2 e per valori di n da 1 a 22. Come mostra Tabella 2 otteniamo
dei buoni risultati e la discrepanza massima tra l’energia calcolata e quella di
riferimento risulta ∼ 10−12EHa.
17
n l j κ E [eV] E0 [eV] ∆E/E0
1 0 0.5 -1 -13.605872851050 -13.605691715850 -3.62 10−4
2 0 0.5 -1 -3.4014795337880 -3.4014229289620 -4.53 10−4
2 1 0.5 1 -3.4014795337880 -3.4014229289620 -4.53 10−4
2 1 1.5 -2 -3.4014342496860 -3.4014229289620 -9.06 10−5
3 0 0.5 -1 -1.5117636500900 -1.5117435239830 -3.62 10−4
3 1 0.5 1 -1.5117636500900 -1.5117435239830 -3.62 10−4
3 1 1.5 -2 -1.5117502325640 -1.5117435239830 -1.21 10−4
3 2 1.5 2 -1.5117502325640 -1.5117435239830 -1.21 10−4
3 2 2.5 -3 -1.5117457601670 -1.5117435239830 -4.03 10−5
4 0 0.5 -1 -0.8503649304807 -0.8503557322406 -2.94 10−4
4 1 0.5 1 -0.8503649304807 -0.8503557322406 -2.94 10−4
4 1 1.5 -2 -0.8503592699680 -0.8503557322406 -1.13 10−4
4 2 1.5 2 -0.8503592699680 -0.8503557322406 -1.13 10−4
4 2 2.5 -3 -0.8503573831737 -0.8503557322406 -5.28 10−5
4 3 2.5 3 -0.8503573831737 -0.8503557322406 -5.28 10−5
4 3 3.5 -4 -0.8503564397823 -0.8503557322406 -2.26 10−5
5 0 0.5 -1 -0.5442325954867 -0.5442276686340 -2.46 10−4
5 1 0.5 1 -0.5442325954867 -0.5442276686340 -2.46 10−4
5 1 1.5 -2 -0.5442296973077 -0.5442276686340 -1.01 10−4
5 2 1.5 2 -0.5442296973077 -0.5442276686340 -1.01 10−4
5 2 2.5 -3 -0.5442287312689 -0.5442276686340 -5.31 10−5
5 3 2.5 3 -0.5442287312689 -0.5442276686340 -5.31 10−5
5 3 3.5 -4 -0.5442282482522 -0.5442276686340 -2.90 10−5
5 4 3.5 4 -0.5442282482522 -0.5442276686340 -2.90 10−5
5 4 4.5 -5 -0.5442279584429 -0.5442276686340 -1.45 10−5
Tabella 1: V (r) = −1/r : Gli autovalori calcolati con il metodo numeri-
co discusso nel testo confrontati con gli autovalori non relativistici E0 di
Eq. (3). L’ultima colonna riporta la differenza relativa tra il valore numerico
relativistico e quello relativo non relativistico.
18
n E [EHa] Eref [EHa] E − Eref [EHa]
1 -0.500006656596500 -0.50000665659536 -1.14 10−12
2 -0.125002080189200 -0.12500208018901 -1.94 10−13
3 -0.055556295176420 -0.05555629517767 1.25 10−12
4 -0.031250338029120 -0.03125033803008 9.57 10−13
5 -0.020000181058520 -0.02000018105900 4.84 10−13
6 -0.013888996749740 -0.01388899675294 3.20 10−12
7 -0.010204150942830 -0.01020415094536 2.53 10−12
8 -0.007812547128859 -0.00781254713001 1.16 10−12
9 -0.006172872986719 -0.00617287298519 -1.53 10−12
10 -0.005000024629192 -0.00500002462650 -2.69 10−12
11 -0.004132250045623 -0.00413225004740 1.78 10−12
12 -0.003472236667825 -0.00347223666540 -2.42 10−12
13 -0.002958591301777 -0.00295859129983 -1.95 10−12
14 -0.002551029591735 -0.00255102959272 9.88 10−13
15 -0.002222229716981 -0.00222222971570 -1.28 10−12
16 -0.001953131195804 -0.00195313119548 -3.26 10−13
17 -0.001730108986653 -0.00173010899016 3.51 10−12
18 -0.001543214251831 -0.00154321425362 1.78 10−12
19 -0.001385045279934 -0.00138504528149 1.56 10−12
20 -0.001250003203452 -0.00125000320622 2.77 10−12
21 -0.001133789620469 -0.00113378962124 7.66 10−13
22 -0.001033060266565 -0.00103306026722 6.53 10−13
Tabella 2: Per i primi 22 stati S (cioe l = 0, j = 1/2, κ = −1), si elenca
l’energia calcolata numericamente nel presente lavoro ed il valore Eref di
energia riportato in letteratura scientifica [8]. L’ultima colonna riporta la
discrepanza, che e compatibile con l’errore numerico del presente calcolo.
19
Figura 2: Autovalori di energia calcolati con potenziale Coulombiano
V = −1/r per n = 5. I dati relativi sono riportati in Tabella 1. Gli split-
ting rispetto a j e la degenerazione rispetto a l in accordo con l’equazione
relativistica all’ordine α2: Eq. (53).
E inoltre possibile verificare che il contributo relativistico ∆E = E−E0 allo
spettro cresce in ottima approssimazione come Z4. Figura 3 mostra il contributo
∆E per il primo livello energetico n = 1, l = 0, j = 1/2 fittato con un andamento
Z4. Come si puo notare abbiamo una buona correlazione tra i valori di ∆E
calcolati e l’andamento Z4 previsto dalla teoria delle perturbazioni al secondo
ordine in α. Lo scostamento dell’ultimo punto dall’equazione (53) e dovuto alle
correzioni di ordine superiore ad α2 che vengono prese in considerazione nella
teoria di Dirac.
20
1 10 100Z
10-5
10-4
10-3
10-2
10-1
100
101
102
103
104
∆E [
EH
a]
Figura 3: Andamento del contributo relativistico ∆E (cerchi) allo sta-
to |nlj 〉 = |1, 0, 1/2 〉. La curva continua rappresenta il risultato dell’e-
quazione relativistica al second’ordine (53). Come si puo notare c’e una
buona correlazione tra i valori calcolati e l’andamento Z4.
4.2 Potenziale Coulombiano a nucleo finito
Consideriamo il potenziale Coulombiano derivante da una distribuzione di carica
nucleare finita. Il potenziale prodotto da una distribuzione sferica di carica di
estensione finita puo essere espresso come
V (r) = −Zeff(r)
r(54)
Zeff(r) = Z1r + Z2r2 + Z3r
3 + ... (55)
Esistono differenti modelli per Zeff(r) reperibili in letteratura. Noi consideriamo
[9]
V (r) =
{− 3Z
2rN
(1 − r2
3r2
N
)se 0 ≤ r ≤ rN
−Zr
se r > rN
(56)
in cui Z1 = 3Z2rN
, Z3 = − Z2r3
N
e tutti gli altri coefficenti Zi sono nulli, che descrive
una distribuzione uniforme di raggio rN . I coefficenti Zi sono fissati richiedendo
21
0,0 2,5×10-5
5,0×10-5
r [a0]
-140000
-120000
-100000
-80000
-60000
-40000
-20000
0
v(r)
[E
Ha]
rN
Figura 4: Confronto tra il potenziale con nucleo puntiforme −Ze2
r (curva
sottile) e il potenziale (56) calcolato per una distribuzione sferica uniforme
di carica di raggio rN (curva spessa). Per r > rN (linea sottile verticale)
i due potenziali coincidono, mentre per r < rN il nucleo finito prevede un
potenziale parabolico.
che V (r) sia continuo con derivata continua in r = rN . Per rN assumo la taglia
del nucleo d’idrogeno (protone) rN = 1.6569 10−5a0 [10]. Figura 4 mostra che
nell’origine questo potenziale raggiunge parabolicamente un valore finito.
La Tabella 3 mostra che il nucleo finito induce una rottura della degenera-
zione rispetto a l. Questo effetto di fatto molto piccolo contribuisce, insieme ad
altri termini piu rilevanti, al Lamb Shift [2].
4.3 Effetti relativistici in presenza di screening
Considero ora un potenziale di tipo Coulombiano che tenga conto dell’attrazione
del nucleo e allo stesso tempo della repulsione degli altri elettroni. Un potenziale
di questo tipo si puo scrivere come (54) in cui
Zeff(r) =Z + r2[α1r + α2r
2 + ...+ αnrn + α0(q + 1)(e
rβ0 − 1)]
1 + r2[β1r + β2r2 + ... + βnrn + α0(er
β0 − 1)](57)
22
n l j κ E [eV]
1 0 0.5 -1 -13.60587284806
2 0 0.5 -1 -3.401479533414
2 1 0.5 1 -3.401479533788
2 1 1.5 -2 -3.401434249686
3 0 0.5 -1 -1.511763649979
3 1 0.5 1 -1.511763650090
3 1 1.5 -2 -1.511750232564
3 2 1.5 2 -1.511750232564
3 2 2.5 -3 -1.511745760167
Tabella 3: In evidenza la rottura di degenerazione rispetto a l causata
dal nucleo finito. Per gli stati 3P3/2 e 3D3/2 lo splitting risulta minore di
10−12 eV e non e quindi evidenziato in tabella.
invece di (55) e dove il numero atomico Z rappresenta la carica nucleare, q la
carica dello ione e i coefficenti αi e βi permettono di interpolare un potenziale
reale. Il termine r2, posto prima delle parentesi quadre, garantisce che il primo
termine correttivo V1(r) al potenziale Coulombiano abbia, per r piccolo, l’an-
damento di r2 , come ci si aspetta da una distribuzione di carica sviluppabile
in serie di potenze in r = 0 e che tenda ad una costante finita nell’origine. In
particolare il termine β0 caratterizza la distanza, in unita di raggi di Bohr, alla
quale la correzione dovuta all’esponenziale risulta significativa. In effetti questa
funzione in un intorno di zero e per distanze molto grandi dall’origine si riduce ai
potenziali Coulombiani −Zr
e − q+1r
rispettivamente. Possiamo cosı imporre che
la soluzione che stiamo cercando soddisfi le medesime condizioni al contorno delle
autofunzioni di un atomo idrogenoide.
Figura 5 mostra l’andamento della carica effettiva definita in Eq. (57) con
α0 = 0.2, β0 = 2.5 e tutti gli altri coefficenti αi e βi nulli e considerando q = 0
e Z = 10. Questi parametri rappresentano una possibile interpolazione del
potenziale dell’atomo di Ne.
Le correzioni relativistiche ad ordine α2, indicate in Eq. (4) si possono
concettualmente dividere in 2 termini:
Hrel = Hs−o −p4
8µ3c2(58)
dove Hs−o tiene conto dell’accoppiamento spin-orbita dell’elettrone e il secondo
termine e la correzione relativistica all’energia cinetica della particella.
23
0 1 2 3 4 5 6 7 8r [a
0]
0
2
4
6
8
10Z
eff(r
)
Figura 5: La carica effettiva Zeff(r) definita dall’equazione (57) “vista” ad
una distanza r dal nucleo per un Atomo di Ne (Z = 10). Si puo notare che
vicino al nucleo (r → 0) l’elettrone non risente della presenza degli elettroni
piu esterni e Zeff(r) → Z, mentre lontano dal nucleo ed entro pochi raggi di
Bohr Zeff(r) → 1 e l’elettrone “sente” il potenziale −1/r tipico dell’atomo
di idrogeno.
Nello spettro dell’Hamiltoniana H e possibile notare uno splitting di energia
tra i livelli j a fissato n e l. Questo e dovuto essenzialmente al termine spin-orbita
dell’Hamiltoniana relativistica. Infatti possiamo scrivere [2]
Hs−o = ξ~S · ~L
~2(59)
dove il parametro ξ di energia di spin-orbita nel caso di potenziale Coulombiano
(52) e
ξ = Z4α2EHa1
n3l(l + 1)(2l + 1)(60)
e l’autovalore dell’operatore ~S · ~L:
⟨l, s, j,mj|
~S · ~L
~2|l, s, j′, m′
j
⟩=j(j + 1) − s(s+ 1) − l(l + 1)
2δjj′δmjm′
j. (61)
Andando quindi a calcolare lo splitting di energia per j = l ± 12
si puo calco-
lare il parametro ξ e studiarne l’andamento rispetto a Z nel caso del potenziale
24
stato n l n3l(l + 1)(2l + 1) Zeff,s−o Ne
2p 2 1 48 10.1177
3p 3 1 162 10.7100
4p 4 1 384 10.6909
5p 5 1 750 7.2277
3d 3 2 810 10.7911
4d 4 2 1920 11.1553
5d 5 2 3750 7.0372
Tabella 4: Importanza del termine n3l(l + 1)(2l + 1) che compare in
Eq. (60). A parita di numero atomico Z, ξ risulta maggiore negli stati
p rispetto agli stati d. Questo significa che gli effetti di spin orbita sono
maggiori negli stati piu vicini al nucleo, dove la velocita dell’elettrone si
avvicina maggiormente a c. La quinta colonna rappresenta il valore di
Zeff,s−o calcolato a partire dall’equazione (62) per l’atomo di Ne.
Coulombiano −Ze2
r(qui ci aspettiamo ξ ∝ Z4) e nel caso del potenziale efficace
Veff di equazioni (54) e (57).
Figura 6 mostra l’andamento di ξ ∝ Z4 aspettato nel caso di potenziale
Coulombiano.
Per un potenziale qualunque V (r), se si pretende di applicare Eq. (60) con
un valore efficace della carica nucleare Zeff,s−o appropriato per ciascun orbitale, e
possibile calcolare Zeff,s−o in termini dello splitting osservato tra gli stati j = l± 12,
per ciascun n:
Zeff,s−o =
((E(l+ 1
2) −E(l− 1
2))
2n3l(l + 1)
α2EHa
) 1
4
(62)
in cui Zeff,s−o e utilizzato al posto di Z in Eq. (60) nel caso sull’elettrone agisca un
potenziale efficace non Coulombiano, ad esempio quello definito dalle equazioni
(54) e (57). Si nota, inoltre, che ξ (60) dipende anche da n3l(l + 1)(2l + 1)
e come si vede in Tabella 4 questo fattore determina l’importanza del termine
ξ per i diversi stati |n, l〉. Diversamente da quanto ci si aspetta, per alcuni
livelli otteniamo Zeff,s−o > Z. Questo e probabilmente dovuto all’utilizzo di una
parametrizzazione (57) del potenziale efficace eccessivamente semplificata.
Figura 7 mostra l’andamento del parametro ξ in funzione del numero atomico
Z nel caso di potenziale efficace di equazioni (54) e (57). Si puo notare che, in
particolare per gli orbitali 2p e 3p maggiormente concentrati vicino al nucleo
25
1 10 100Z
10-8
10-7
10-6
10-5
10-4
10-3
10-2
10-1
100
101
102
ξ(Ζ) 2p
3p3d4p4d4f5p5d5f5g
Figura 6: Andamento del parametro calcolato di spin orbita ξ (punti)
in funzione del numero atomico Z per gli stati 2p, 3p, 3d, 4p, 4d, 4f , 5p,
5d, 5f e 5g nel caso di potenziale Coulombiano (52). Le curve continue o
tratteggiate rappresentano il risultato dell’espressione (60). Come si nota
c’e un ottimo accordo tra i valori calcolati e l’andamento della formula
perturbativa di ordine α2. In particolare per Z = 1, il nostro metodo di
calcolo ci conduce a ξ2p = 0.030009 meV, in buon accordo con la predizione
di Eq. (60) ξ2p,Eq. = 0.030188 meV (∆ ≃ 0.179 µeV).
26
1 10 100Z
10-6
10-4
10-2
100
102
ξ(Ζ)
2p3p3d4p4d5p5d
Figura 7: Andamento del parametro di spin orbita ξ (punti) in funzione
del numero atomico Z per gli stati 2p, 3p, 3d, 4p, 4d, 5p e 5d nel caso di
potenziale efficace definito dalle equazioni (54) e (57).
27
(bassi valori di r), l’andamento e simile a quello di Figura 6. Infatti per r → 0
si puo notare dall’equazione (57) e dalla Figura 5 che Zeff(r) → Z e quindi i due
potenziali sono simili e cosı i loro spettri.
5 Discussioni e Conclusioni
Nel presente lavoro di tesi abbiamo implementato un programma di calcolo rela-
tivamente semplice per risolvere numericamente l’equazione di Dirac radiale per
un singolo elettrone in moto in un potenziale centrale V (r) qualunque. Abbia-
mo testato il codice sul caso ben noto del potenziale Coulombiano prodotto da
un nucleo puntiforme (52). Abbiamo poi considerato una distribuzione di carica
nucleare finita (56), nonche una semplice forma di potenziale efficace per il Ne,
Eq. (54) e (57). Abbiamo potuto osservare gli spettri nei diversi casi, individuan-
do vari effetti interessanti negli spettri, soffermandoci in particolare sugli effetti
relativistici.
L’equazione di Dirac radiale deve essere risolta imponendo le condizioni
al contorno nell’intervallo [r1, r2) = [0,∞). Dovendo pero risolvere numerica-
mente il problema, e stato necessario integrare l’equazione in uno spazio finito ed
evitare possibili singolarita nell’origine. Di fatto abbiamo integrato l’equazione
tra r1 ∼ 10−12 a0 e r2 ∼ 50 a0. Per gli atomi piu pesanti, al crescere di Z, e
risultato necessario ridurre i valori sia di r1 che di r2: di fatto, per Z grandi,
le funzioni d’onda si concentrano vicino all’origine e vanno rapidamente a zero
entro pochi raggi di Bohr. risulta quindi importante integrare in modo raffinato
nelle zone vicino all’origine, andando a ridurre il passo di integrazione, ove neces-
sario. Abbiamo notato che, per determinare in maniera efficace gli autovalori, le
condizioni al contorno sull’autofunzione gκ contano relativamente poco. In par-
ticolare, nel punto x2 in cui abbiamo imposto che la funzione gκ si annulli, in
realta, per quanto concerne il calcolo degli autovalori, avremmo potuto supporre
semplicemente che gκ(x2) abbia un valore finito.
Il metodo di calcolo utilizzato si e dimostrato soprendentemente efficace.
Abbiamo riscontrato una buona compatibilita tra i risultati da noi ottenuti e
quelli reperibili in letteratura scientifica.
28
Riferimenti bibliografici
[1] J. J. Sakurai, Meccanica Quantistica Moderna (Zanichelli, Bologna, 2003).
[2] N. Manini, Introduction to the Physics of Matter (Cusl, Milano, 2008).
[3] P. A. Dirac, Proc. R. Soc. London A 117, 610-624 (1928).
[4] O. Certık, Numerical solution of the radial Dirac equation in pseudopotential
construction (Bachelor Thesis, Charles University in Prague, 2007),
http://ondrej.certik.cz/site_media/cookbook/thesis.ps
[5] F. Salvat and R. Mayol, Comp. Phys. Comun. 62, 65-79 (1991).
[6] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Nu-
merical Recipies in C++. The Art of Scientific Computing. Second Edition
(Cambridge University Press, 2002).
[7] G. Naldi, L. Pareschi e G. Russo, Introduzione al Calcolo Scientifico, metodi
e applicazioni con Matlab (McGraw-Hill, Milano, 2003).
[8] M. I. Bhatti and W. F. Perger, J. Phys. B: At. Mol. Opt. Phys. 39, 553-558
(2006).
[9] I. P. Grant, J. Phys. B: At. Mol. Opt. Phys. 42, 1-8 (2009).
[10] Fundamental Physical Constant from NIST:
http://physics.nist.gov/cuu/Constants/index.html
29
Ringraziamenti
Non avrei mai creduto che ce l’avrei fatta. Se oggi sono qui, felice per aver
raggiunto questa tappa importante della mia vita, certamente devo ringraziare
tutti coloro che mi hanno sostenuto e aiutato. Innanzitutto un doveroso grazie
a Nicola Manini per la sua pazienza e disponibilita in questi tre mesi di lavoro,
che mi hanno aiutato ad amare ancora di piu questa scienza: la Fisica. Un grazie
al professor Favale che ha provato a spiegarmela ai tempi del Liceo e che mi ha
trasmesso la sua grande passione per i ragazzi e per l’insegnamento.
Un grazie alle persone con le quali ho condiviso questi tre anni di studio
e grazie alle quali sono riuscito a passare anche qualche esame: a J. e al Sanze
per la loro infinita disponibilita a spiegarmi le cose, a Marta ed Emy per la loro
capacita di sorridere sempre e comunque, a Giulia, Alessio, Anna e Paolo.
Un grazie speciale e per la mia meravigliosa famiglia: a mio papa e a mia
mamma, genitori stupendi, a Marta, Annalisa e Gabriele. A nonna Diletta, a
Davide, a Giada e alla piccola Rebecca, bella come il sole! Questo grazie e per
dirvi (anche se lo faccio raramente) che vi voglio proprio un sacco di bene.
Un grazie va agli Amici del Sidamo ed in modo particolare al meraviglioso
gruppo di “Cernusco on the River”, veri compagni di viaggio, perche mi aiutano
a crescere nel servizio e nel “faticare di brutto” per i ragazzi e per i poveri. E un
grazie agli amici di Abobo che porto sempre nel cuore e a tutti gli amici volontari,
che mi richiamano ad essere una persona semplice e a mettere al centro della mia
vita cio che conta davvero: “Amare, distribuire speranza, fiducia e gioia”.
Un grazie particolare a Don Giuliano e a Don Fiorenzo, preziose guide nel
cammino per arrivare al Signore.
Un ringraziamento particolare ai miei Amici, quelli con la “A” maiuscola
con i quali condivido sogni e ideali e vorrei continuare a farlo. Grazie alla Manu
sempre umile e silenziosa, alla Ceci per il suo grande desiderio di voler bene, alla
Pannalis per la sua enorme attenzione agli altri, all’Ari per la sua folle pazzia e il
suo amore per i giovani, a Joel per il suo entusiasmo travolgente nel fare le cose,
a Paolo per l’allegria con la quale sa vivere ogni giorno, alla Rossi per il desiderio
di condividere la sua voglia di servire, al Giova per il suo coraggio di amare in
modo radicale e fino in fondo. Grazie perche siete persone cosı speciali per la mia
vita.
30
Un grazie ad Anna per l’enorme amore che mette in tutto cio che fa, per la
sua capacita di amarmi per quello che sono, per avermi insegnato a sognare e per
il suo desiderio di volare, e volare alto, nella vita.
“Vivi la vita, sogna il tuo sogno, chiamalo Amore”
...un ringraziamento al Signore perche ha messo proprio voi sul mio cammino
e perche mi permette di sognare in grande...
31