Teoria dei condensati di Bose-Einstein

Appunti del corso di Fenomeni Quantistici Macroscopici

Lezioni tenute dal Dott. M. Modugno

A.A. 2008-2009 / 2009-2010

Ringrazio Giacomo Cappellini (A.A. 2008/2009) che ha provveduto allaprima stesura di questi appunti, e Vladislav Gavryusev (A.A. 2009/2010) persuggerimenti e correzioni.

Indice

Introduzione 1

1 Introduzione alla teoria dei campi 3Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 Il potenziale interatomico 112.1 Il metodo dello pseudo-potenziale . . . . . . . . . . . . . . . . 12

2.1.1 Modello a due sfere rigide . . . . . . . . . . . . . . . . 122.1.2 Modello con potenziale interatomico . . . . . . . . . . 142.1.3 Problema a N corpi . . . . . . . . . . . . . . . . . . . . 15

Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3 La teoria di Gross-Pitaevskii 193.1 L’equazione di Gross-Pitaevskii . . . . . . . . . . . . . . . . . 203.2 L’equazione di Gross-Pitaevskii stazionaria . . . . . . . . . . . 21

3.2.1 Il potenziale armonico . . . . . . . . . . . . . . . . . . 233.2.2 Il limite Thomas-Fermi . . . . . . . . . . . . . . . . . . 24

3.3 Formulazione idrodinamica . . . . . . . . . . . . . . . . . . . . 26Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

4 Equazioni di scaling 294.1 Scaling delle equazioni idrodinamiche . . . . . . . . . . . . . . 294.2 Scaling dell’equazione di Gross-Pitaevskii . . . . . . . . . . . . 32

4.2.1 Scaling Thomas-Fermi . . . . . . . . . . . . . . . . . . 344.2.2 Scaling non-interagente . . . . . . . . . . . . . . . . . . 34

4.3 Soluzione delle equazioni di scaling . . . . . . . . . . . . . . . 354.3.1 Espansione libera . . . . . . . . . . . . . . . . . . . . . 354.3.2 Piccole oscillazioni . . . . . . . . . . . . . . . . . . . . 40

4.4 Scaling di un gas termico non interagente . . . . . . . . . . . . 424.4.1 Distribuzione di equilibrio . . . . . . . . . . . . . . . . 434.4.2 Espansione libera . . . . . . . . . . . . . . . . . . . . . 44

4.5 Distribuzione dei momenti . . . . . . . . . . . . . . . . . . . . 46Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

i

5 Eccitazioni collettive 495.1 Eccitazioni di Bogoliubov . . . . . . . . . . . . . . . . . . . . 49

5.1.1 Eccitazioni di un condensato uniforme . . . . . . . . . 525.1.2 Criterio di Landau per la superfluidità . . . . . . . . . 55

5.2 Vortici quantizzati . . . . . . . . . . . . . . . . . . . . . . . . 56Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

6 Condensati in reticoli ottici 596.1 L’hamiltoniana di Bose-Hubbard . . . . . . . . . . . . . . . . 59

6.1.1 Approssimazione gaussiana . . . . . . . . . . . . . . . . 636.1.2 Calcolo dell’energia di interazione U . . . . . . . . . . . 646.1.3 Calcolo dell’energia di tunneling J . . . . . . . . . . . . 646.1.4 Approssimazione di Gutzwiller e diagramma di fase . . 66

6.2 Instabilità energetica e dinamica . . . . . . . . . . . . . . . . . 686.2.1 Instabilità energetica . . . . . . . . . . . . . . . . . . . 706.2.2 Instabilità dinamica . . . . . . . . . . . . . . . . . . . . 73

Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

ii

Introduzione

I condensati di Bose-Einstein, prodotti in laboratorio utilizzando atomi raf-freddati alle temperature del nanokelvin, rappresentano uno strumento al-tamente versatile per studiare e riprodurre su scala macroscopica un’am-pia gamma di fenomeni quantistici e per simulare il comportamento di altrisistemi fisici.

In primo luogo, infatti, i condensati possono essere considerati come ondecoerenti di materia, formata da atomi che occupano tutti lo stesso stato, ea causa della loro natura intrinsecamente quantistica e delle dimensioni finoalla scala del millimetro, rappresentano un sistema ideale per studiare effettiquantistici su scala macroscopica. Un esempio di questo è il fatto che ladistribuzione di probabilità quantistica diventa di fatto la distribuzione didensità degli atomi del condensato, che può essere misurata tramite immaginifotografiche digitali.

Un’altro aspetto importante riguarda il fatto che manipolando oppor-tunamente i potenziali a cui sono soggetti gli atomi è possibile riprodurresituazioni in cui l’hamiltoniana che descrive il sistema è la stessa di un altrosistema fisico, di fatto quindi simulando il comportamento di quest’ultimo(un esempio ne sono i condensati in reticoli ottici, che posso simulare - percerti aspetti - il comportamento degli elettroni in un cristallo). Questo aspet-to realizza l’idea avuta da Feynman agli inizi degli anni ’80 di simulatorequantistico, e rappresenta una delle attuali frontiere della fisica degli atomiultrafreddi.

In questi appunti verranno introdotti alcuni concetti generali che stannoalla base della descrizione teorica dei condensati di Bose-Einstein.

1

2

Capitolo 1

Introduzione alla teoria deicampi

In questo capitolo mostriamo l’equivalenza tra Meccanica Quantistica perun sistema di N particelle e Teoria dei Campi. Quest’ultima formulazionepermette di descrivere il sistema in modo compatto ed efficiente tramiteoperatori di campo che dipendono da un solo set di coordinate spaziali ~rinvece che dalle coordinate {~r1, ~r2 . . . , ~rN} di tutte le particelle.

Consideriamo un gas di N particelle di massa m interagenti con intera-zioni a due corpi, in presenza di un potenziale esterno V (~r). Questo sistemaè descritto dall’hamiltoniana

H = − ~2

2m

∑i

∇2i +∑i

V (~ri) +∑i

Le autofunzioni ψn e lo spettro di autovalori En dell’hamiltoniana in (1.1)si ottengono risolvendo l’equazione di Schrödinger stazionaria

Hψn(~r1, . . . , ~rN) = Enψn(~r1, . . . , ~rN) (1.3)

dove la funzione d’onda ψn è normalizzata nel seguente modo∫d3Nr |ψ|2 = 1

La risoluzione dell’equazione diventa impraticabile quando N è un numerogrande. Conviene quindi utilizzare la teoria dei campi. Faremo di seguito unabreve introduzione assiomatica alla teoria dei campi, seguendo la trattazionein [1].

Nella teoria dei campi la funzione d’onda - che dipende da 3N coordinatespaziali - viene sostituita da un operatore di campo ψ̂ che dipende solo da 3coordinate spaziali

ψ(~r1, . . . , ~rN) → ψ̂(~r). (1.4)In questo modo si passa ad una rappresentazione che tiene conto della totalitàdel sistema e non delle singole particelle che lo compongono.

In questo capitolo mostreremo che questa formulazione è del tutto analogaalla meccanica quantistica. Vedremo poi che nel caso dei condensati di Bose-Einstein è possibile - nell’approssimazione di campo medio - ritornare ad unadescrizione tipo Schrödinger in termini di una funzione d’onda (non più uncampo) che dipende da un solo set di tre coordinate spaziali, e che riveste ilruolo di funzione d’onda macroscopica del condensato.

Gli operatori di campo sono definiti tramite le seguenti regole di commu-tazione (nel caso di bosoni):[

ψ̂(~r), ψ̂†(~r ′)]

= δ(~r − ~r ′) (1.5)[ψ̂(~r), ψ̂(~r ′)

]= 0[

ψ̂†(~r), ψ̂†(~r ′)]

= 0

Nel caso di fermioni i commutatori sono sostituiti dagli anticommutatori{ψ̂(~r), ψ̂†(~r ′)

}= δ(~r − ~r ′) (1.6){

ψ̂(~r), ψ̂(~r ′)}

= 0{ψ̂†(~r), ψ̂†(~r ′)

}= 0

Analogamente, l’hamiltoniana è ottenuta tramite operatori costruiti intermini dell’hamiltoniana di partenza e degli operatori di campo, nel seguentemodo

Ĥ = K̂ + V̂ + Ω̂ (1.7)

4

con

K̂ = − ~2

2m

∫d~r ψ̂†(~r)∇2ψ̂(~r) (1.8)

V̂ =

∫d~r ψ̂†(~r)V (~r)ψ̂(~r) (1.9)

Ω̂ =1

2

∫d~r1

∫d~r2 ψ̂

†(~r1)ψ̂†(~r2)v(~r1, ~r2)ψ̂(~r1)ψ̂(~r2) (1.10)

dove il fattore 1/2 nella terza equazione serve per contare correttamente leinterazioni a due corpi.

Introduciamo ora l’operatore numero N̂ , definito come

N̂ =

∫d~r ψ̂†(~r)ψ̂(~r) (1.11)

che, come vedremo tra poco, è associato al numero totale di particelle delsistema. Si dimostra facilmente che l’operatore commuta con l’hamiltoniana(questa proprietà vale sia per i bosoni che per i fermioni)[

Ĥ, N̂]

= 0 (1.12)

ed è quindi possibile scegliere una base di autostati simultanei per Ĥ e N̂ .Indicheremo tali autostati con |ψEN〉. Questo autostati sono tali che

Ĥ |ψEN〉 = E |ψEN〉 (1.13)N̂ |ψEN〉 = N |ψEN〉 (1.14)

un’altra proprietà dei vettori |ψEN〉 è che

〈ψEN |ψEN〉 = 1

un particolare stato della base |ψEN〉 è lo stato |ψ00〉, che indicheremo con|0〉, e che prende il nome di stato di vuoto. Questo è un autostato deglioperatori Ĥ e N̂ con autovalore nullo

Ĥ |0〉 = 0 N̂ |0〉 = 0.

Conoscendo la forma dell’operatore numero e tenendo conto del commutatore(1.5) è possibile dimostrare che[

ψ̂(~r), N̂]

= ψ̂(~r) (1.15)[ψ̂†(~r), N̂

]= −ψ̂†(~r) (1.16)

5

Possiamo vedere come questi commutatori presentino un’analogia con i com-mutatori tra gli operatori di creazione e distruzione ↠e â e l’operatore dioccupazione n̂ = â†â nel caso dell’oscillatore armonico. Procedendo come intale caso, possiamo mostrare come gli operatori ψ̂† e ψ̂ operano sugli stati|ψEN〉. Ad esempio, usando la (1.15) si ha che(

ψ̂N̂ − N̂ψ̂)|ψEN〉 = ψ̂ |ψEN〉

ed analogamente per ψ̂† (con il segno meno), da cui si ottiene:

N̂ψ̂ |ψEN〉 = (N − 1)ψ̂ |ψEN〉N̂ψ̂† |ψEN〉 = (N − 1)ψ̂† |ψEN〉

Procedendo come nel caso dell’oscillatore armonico, dalle equazioni prece-denti è possibile mostrare che gli autovalori dell’operatore N̂ sono numeriinteri e che

ψ̂ |ψEN〉 =√N |ψEN−1〉 (1.17)

ψ̂† |ψEN〉 =√N + 1 |ψEN+1〉 (1.18)

Queste equazioni mostrano che gli operatori ψ̂† e ψ̂ sono degli operatoririspettivamente di creazione e distruzione, analogamente agli operatori ↠eâ nel caso dell’oscillatore armonico. In particolare ψ̂† e ψ̂ sono operatoridi creazione e distruzione di particelle, fanno cioè passare da stati con nparticelle a stati con N+1 o N−1 particelle (nell’oscillatore armonico creanoe distruggono un quanto di energia facendo passare il sistema - una singolaparticella che non viene né creata né distrutta - da un livello energetico aduno adiacente).

A questo punto possiamo mostrare come è possibile costruire tramitegli operatori di campo la funzione d’onda ad N particelle della meccanicaquantistica. Notiamo intanto che, per la (1.17), vale

〈ψE′N ′| ψ̂(~r1), . . . , ψ̂(~rN) |ψEN〉 ∝ 〈ψE′N ′|0〉

e supponiamo che lo stato di vuoto sia unico

〈ψE′N ′|0〉 6= 0 ⇔ |ψE′N ′〉 = |0〉 .

Possiamo quindi costruire una funzione del tipo

ψEN(~r1, . . . , ~rN) =1√N !〈0| ψ̂(~r1) . . . ψ̂(~rN) |ψEN〉 (1.19)

Nel seguito dimostreremo che questa funzione ha tutte le caratteristiche dellafunzione d’onda, cioè è correttamente normalizzata e soddisfa l’equazione di

6

Schrödinger stazionaria (1.3). Verifichiamo innanzitutto che la funzione ènormalizzata∫d3Nr ψ∗EN(~r1, . . . , ~rn)ψEN(~r1, . . . , ~rN) =

=1

N !

∫d3Nr 〈ψEN | ψ̂†(~rN) . . . ψ̂†(~r1) |0〉 〈0| ψ̂(~r1) . . . ψ̂(~rN) |ψEN〉 =

=1

N !

∑E′,N ′

∫d3Nr 〈ψEN | ψ̂†(~rN) . . . ψ̂†(~r1) |ψE′N ′〉 〈ψE′N ′| ψ̂(~r1) . . . ψ̂(~rN) |ψEN〉 =

=1

N !

∫d3Nr 〈ψEN | ψ̂†(~rN) . . . ψ̂†(~r1)ψ̂(~r1) . . . ψ̂(~rN) |ψEN〉 =

=1

N !

∫d3(N−1)r 〈ψEN | ψ̂†(~rN) . . . ψ̂†(~r2)N̂ψ̂(~r2) . . . ψ̂(~rN) |ψEN〉

Nei passaggi precedenti sono stati prima aggiunti i contributi delle varie〈ψEN | ψ̂(~r1), . . . , ψ̂(~rN) |ψE′N ′〉 che sono nulli per |ψE′N ′〉 6= |0〉, e poi sfruttatala proprietà di completezza:∑

E′N ′

|ψE′N ′〉 〈ψE′N ′| = 1

Utilizzando poi la regola di commutazione (1.15) si ha

N̂ψ̂ = ψ̂(N̂ − 1)

e riprendendo dall’ultima equazione del blocco si ha che la catena di ugua-glianze continua con

1

N !

∫d3(N−1)r 〈ψEN | ψ̂†(~rN) . . . ψ̂†(~r2)ψ̂(~r2)(N̂ − 1) . . . ψ̂(~rN) |ψEN〉 =

1

N !

∫d3(N−2)r 〈ψEN | ψ̂†(~rN) . . . ψ̂†(~r3)ψ̂(~r3)N̂(N̂ − 1) . . . ψ̂(~rN) |ψEN〉 ;

iterando il procedimento di semplificazione utilizzando i commutatori si ot-tiene

1

N !〈ψEN | N̂ · (N̂ − 1) · · · · 1 |ψEN〉 = 1.

Vediamo ora di dimostrare che la funzione ψEN(~r1, . . . , ~rN) cos̀ı com’èstata costruita soddisfa il seguente teorema che completa l’analogia tra campoquantizzato e sistema ad N corpi (per brevità di notazione definiamo Vi =V (~ri), vij = v(~ri, ~rj))

Teorema 1. Sia ψEN(~r1, . . . , ~rN) =1√N !〈0| ψ̂(~r1) . . . ψ̂(~rN) |ψEN〉. Questa

soddisfa l’equazione di Schrödinger(− ~

2

2m

∑i

∇2i +∑i

Vi +∑i

Dimostrazione. Dalla (1.13) e dalla definizione di ψEN(~r1, . . . , ~rN) si ha che

1√N !〈0| ψ̂(~r1) . . . ψ̂(~rN)Ĥ |ψEN〉 = EψEN(~r1, . . . , ~rN) (1.20)

alla precedente possiamo sottrarre il termine 〈0| Ĥψ̂(~r1) . . . ψ̂(~rN) |ψEN〉, cheè nullo essendo Ĥ |0〉 = 0. La (1.20), può essere quindi riscritta come

1√N !〈0|[ψ̂(~r1) . . . ψ̂(~rN), Ĥ

]|ψEN〉 = EψEN(~r1, . . . , ~rN) (1.21)

Il commutatore dell’equazione precedente può essere calcolato usando a ca-tena la regola di commutazione

[ÂB̂, Ĉ] = Â[B̂, Ĉ] + [Â, Ĉ]B̂ (1.22)

con cui si ottiene

1√N !

N∑j=1

〈0| ψ̂(~r1) . . .[ψ̂(~rj), Ĥ

]. . . ψ̂(~rN) |ψEN〉 (1.23)

Calcoliamo il commutatore [ψ̂(~rj), Ĥ]. Dato che Ĥ = K̂ + V̂ + Ω̂ calcoliamoil commutatore per ciascuno dei tre termini.

Per quanto riguarda il termine di energia cinetica si ha

[ψ̂(~rj), K̂] = −~2

2m

∫d~r[ψ̂(~rj), ψ̂

†(~r)∇2ψ̂(~r)]

=

= − ~2

2m

∫d~r[ψ̂(~rj), ψ̂

†(~r)]∇2ψ̂(~r) =

= − ~2

2m

∫d~r δ(~r − ~rj)∇2ψ̂(~r) =

= − ~2

2m∇2j ψ̂(~rj) (1.24)

Analogamente, il termine di energia potenziale è semplicemente dato da

[ψ̂(~rj), V̂ ] = V (~rj)ψ̂(~rj) (1.25)

Per il commutatore con il termine di energia di interazione abbiamo invece

[ψ̂(~rj), Ω̂] =1

2

∫d~r1

∫d~r2

[ψ̂(~rj), ψ̂

†(~r1)ψ̂†(~r2)v12ψ̂(~r1)ψ̂(~r2)

]=

=1

2

∫d~r1

∫d~r2

[ψ̂(~rj), ψ̂

†(~r1)ψ̂†(~r2)

]v12ψ̂(~r1)ψ̂(~r2) =

=1

2

∫d~r1ψ̂

†(~r1)v1jψ̂(~r1)ψ̂(~rj) +1

2

∫d~r2ψ̂

†(~r2)vj2ψ̂(~r2)ψ̂(~rj) =

=

∫d~rψ̂†(~r)v(~r, ~rj)ψ̂(~r)ψ̂(~rj) (1.26)

8

dove per spezzare i due integrali si è usato la proprietà (1.22), la formaesplicita (1.5) del commutatore, ed infine integrato. Nel passaggio successivosi è sfruttata la simmetria dell’interazione a due corpi, vij = vji.

Calcoliamo adesso il contributo del primo commutatore alla (1.23)

1

N !

∑j

(− ~

2

2m∇2j)〈0| ψ̂(~r1) . . . ψ̂(~rj) . . . ψ̂(~rN) |ψEN〉 =

=∑j

(− ~

2

2m∇2j)ψEN (1.27)

Il secondo termine è banalmente∑j

V (~rj)ψEN (1.28)

Infine, per il termine di interazione tra particelle si ha

1

N !

∑j

〈0| ψ̂(~r1) . . .∫d~rψ̂†(~r)v(~r, ~rj)ψ̂(~r)ψ̂(~rj) . . . ψ̂(~rN) |ψEN〉 ; (1.29)

questa espressione può essere semplificata considerando l’operatore

X̂(~rj) =

∫d~rψ̂†(~r)v(~r, ~rj)ψ̂(~r) (1.30)

che gode delle seguenti proprietà

[ψ̂(~ri), X̂(~rj)] = vijψ(~ri) ; X̂ |0〉 = 0

Dalla prima si ha che

ψ̂(~r1) . . . ψ̂(~rj−1)X̂(~rj)ψ̂(~rj) . . . ψ̂(~rN) =

=

(X̂(~rj) +

j−1∑i=1

vij

)ψ(~r1) . . . ψ(~rN); (1.31)

inserendo questa espressione nella (1.29) ed usando la seconda proprietà, siottiene

N∑j=1

j−1∑i=1

vij1

N !〈0| ψ̂(~r1) . . . ψ̂(~rN) |ψEN〉 =

∑i

Bibliografia

[1] Kerson Huang, Statistical Mechanics, second edition, Wiley (1987),Appendice A.

10

Capitolo 2

Il potenziale interatomico

In questa sezione vedremo come è possibile scrivere in forma approssimata ilpotenziale di interazione a due corpi nel caso di un gas diluito. Questa con-dizione si verifica quando sia la distanza media tra particelle v1/3 ≡ (V/N)1/3che la lunghezza d’onda termica λT sono entrambe molto maggiori del raggiod’azione a del potenziale. La prima ci dice che la probabilità di trovare dueparticelle all’interno del raggio di interazione di qualunque altra è piccola, eche è lecito trascurare interazioni a più corpi. Questa è una condizione chevale anche classicamente. La seconda è invece strettamente legata alle pro-prietà quantistiche: poiché una particella non può localizzarsi su lunghezzepiù piccole della sua lunghezza d’onda, essa risulta quindi estesa su una scalache è molto più grande del raggio di interazione del potenziale ed i dettaglidel potenziale diventano poco importanti. In altre parole, la particella vedesolo un effetto medio del potenziale.

Con queste ipotesi vedremo che si può semplificare la forma del potenzialedi interazione ponendo

v(~ri, ~rj) =4π~2am

δ(~ri − ~rj) ≡ gδ(~ri − ~rj) (2.1)

che permette di riscrivere il termine dell’hamiltoniana di campo come

Ω̂ =g

2

∫d~rψ̂†(~r)ψ̂†(~r)ψ̂(~r)ψ̂(~r) (2.2)

Vediamo ora in dettaglio come si procede per ottenere questa espressione,iniziando dal caso semplificato di un sistema a due particelle, ed estendendopoi i risultati al caso di N particelle.

11

2.1 Il metodo dello pseudo-potenziale

Iniziamo consideriamo un sistema semplice composto da due particelle dimassa m, descritto dall’hamiltoniana

Ĥ =p̂212m

+p̂222m

+ V (|~r1 − ~r2|) (2.3)

se introduciamo le variabili del centro di massa~R =

~r1 + ~r22

(2.4)

~P = ~p1 + ~p2 → −i~∂

∂ ~R(2.5)

e da quelle relative ~r = ~r1 − ~r2 (2.6)~p = ~p1 − ~p22

→ −i~ ∂∂~r

(2.7)

L’hamiltoniana può essere riscritta come

Ĥ =P̂ 2

2M+p̂2

2µ+ V (~r) (2.8)

e si può separare il moto del centro di massa da quello relativo fattorizzandola funzione d’onda in una parte di onda libera del centro di massa a unaparte relativa

ψ(~r1, ~r2) = ei

~P ·~R~ ψ(~r) (2.9)

con la funzione d’onda relativa che soddisfa l’equazione:[p̂2

2µ+ V (~r)

]ψ(~r) = Eψ(~r) (2.10)

2.1.1 Modello a due sfere rigide

Supponiamo ora che il potenziale di interazione V (~r) sia di tipo sfera rigidadi raggio a

V (~r) =

{0 per r > a (2.11)

∞ per r ≤ a ; (2.12)

la funzione d’onda è quindi identicamente nulla per r ≤ a, mentre per r > al’equazione di Schrödinger si riduce a quella di particella libera con energiaE = ~2k2/2µ: {

(∇2 + k2)ψ(~r) = 0 per r > a (2.13)ψ(~r) = 0 per r ≤ a . (2.14)

12

Vogliamo adesso trovare il modo di eliminare le condizioni al contorno edavere una sola equazione che però abbia soluzione coincidente con quella delproblema originario per r > a. Utilizzeremo il metodo degli pseudopotenzialiche consiste nel costruire un potenziale con elementi fittizi che simulano lecondizioni al contorno consentendoci di eliminarle dal problema.

Iniziamo considerando soluzioni a simmetria sferica, nel limite k → 0. Ilsistema precedente diventa

1

r2d

dr

(r2dψ

dr

)= 0 per r > a (2.15)

ψ(r) = 0 per r ≤ a (2.16)

Il sistema può essere risolto esattamente e si ottiene:

ψ(r) =

{χ(1− a

r

)per r > a (2.17)

0 per r ≤ a (2.18)

dove χ è un coefficiente che dipende dalle condizioni al contorno.Definiamo ora un funzione estesa ψex(r) (per r 6= 0), che soddisfa l’equa-

zione:

(∇2 + k2)ψex(r) = 0 (2.19)

con la condizione al contorno ψex(a) = 0. Si ha che per k → 0:

ψex(r) −−→r→0

χ(1− a

r

)(2.20)

dove

χ =∂

∂r(rψex)

∣∣∣∣r=0

(2.21)

Generalizziamo ora l’equazione (2.19) al fine di includere anche il punto r = 0.Questo si può fare una volta determinato il comportamento di (∇2 + k2)ψexvicino a r = 0. Utilizzando la (2.20) e considerando che siamo nell’ipotesi dik → 0 si ha:

∇2ψex(r) −−→r→0

−4πaδ(~r)χ (2.22)

in cui si è usato la nota proprietà ∇2(1/r) = −4πδ(~r). Ricordando ora la(2.21) si ottiene:

∇2ψex(r) −−→r→0

4πaδ(~r)∂

∂r(rψex) (2.23)

In sostanza, per k → 0 la funzione ψex(r) soddisfa ovunque l’equazione:

(∇2 + k2)ψex(r) = 4πaδ(~r)∂

∂r(rψex) (2.24)

13

Si è quindi introdotto uno pseudopotenziale rappresentato dall’operatoreδ(~r)(∂/∂r)r. Per piccoli k e per r > a la ψex(r) soddisfa la stessa equazionele stesse condizioni a contorno di ψ(r), per cui ψex(r) = ψ(r) per r > a.

Tuttavia questo risultato è stato ottenuto nell’approssimazione di piccolik e di onda s. Per ottenere una soluzione più generale è necessario espanderein onde parziali non sfericamente simmetriche e considerare k arbitrari. Sipuò dimostrare che l’esatto pseudopotenziale di onda s è dato da:

− 4πk cot η0

δ(~r)∂

∂rr (2.25)

dove η0 è lo sfasamento in onda s per il potenziale di sfera rigida:

− 1k cot η0

=tan ka

k= a

[1 +

1

3(ka)2 + · · ·

](2.26)

Sarebbe inoltre necessario aggiungere serie infinita di altri pseudopotenzialiche rappresentano gli effetti di diffusione in onde p, d etc., proporzionali aa2l+1.

Da questo risultati si può vedere che la (2.24) è corretta fino all’ordinea2, cioè se la funzione d’onda ψ(r) e l’autovalore k sono sviluppati in serie dipotenze di a allora la (2.24) da i coefficienti di a a a2 corretti.

2.1.2 Modello con potenziale interatomico

Consideriamo ora il caso di due particelle interagenti per mezzo di un genericopotenziale interatomico v(r) a raggio d’azione finito che non ammetta statilegati. Facendo nuovamente la separazione tra variabili del centro di massae relativa si avrà che il moto relativo è governato dall’equazione:[

−~2∇2

2µ+ v(r)

]ψ(~r) = Eψ(~r) (2.27)

ponendo ~2k2 = 2µE si ottiene:

~2

2µ(∇2 + k2)ψ(~r) = v(r)ψ(~r) (2.28)

con delle assegnate condizioni a contorno per r → ∞. A basse energie èimportante soltanto la diffusione delle onde s, pertanto faremo di nuovo leapprossimazioni di simmetria sferica e k → 0. L’equazione precedente siriduce a:

u′′(r) + k2u(r) =2µ

~2v(r)u(r) (2.29)

dove:u(r) ≡ rψ(r) (2.30)

14

Poiché si è assunto che v(r) abbia un raggio d’azione finito e non abbia statilegati, nel limite r →∞ u(r) tende ad essere una funzione sinusoidale

u(r) −−−→r→∞

u∞(r) (2.31)

conu∞(r) ≡ rψ∞(r) ∝ (sin kr + tan η0 cos kr) (2.32)

dove η0 è lo sfasamento in onda s. Inoltre si ha che per k → 0:

ψ∞(r) ∝1

r(sin kr + tan η0 cos kr) −−→

r→0c

(1 +

tan η0kr

)(2.33)

In generale η0 è funzione di k; per piccoli k si ha il seguente sviluppo:

k cot η0 = −1

a+

1

2k2r0 + · · · (2.34)

dove a è la lunghezza di diffusione e r0 il raggio di interazione efficace.Sostituendo questo sviluppo nella (2.33) si ottiene

ψ∞(r) −−→r→0

1 +tan η0kr

= 1− ar

(2.35)

da cui si ottienerψ∞(r) −−→

r→0r − a (2.36)

A questo punto è facile capire il significato della lunghezza di diffusione a:come si può vedere dalla Fig. 2.1, questa è l’intercetta della funzione d’ondaasintotica rψ∞(r) sull’ascissa. Nel caso del potenziale della sfera rigida lalunghezza di diffusione è il diametro della sfera. In generale a può esseresia positivo che negativo, a seconda che il potenziale sia prevalentementerepulsivo o attrattivo.

2.1.3 Problema a N corpi

Dopo aver introdotto gli pseudopotenziali nel problema a due corpi possiamoora discutere la generalizzazione al problema a N corpi.

Consideriamo innanzitutto il problema a N corpi interagenti come sfererigide, e definiamo rij = |~ri − ~rj|. L’equazione di Schrödinger del sistema è−

~2

2m(∇21 + · · ·+∇2N)Ψ = EΨ se rij > a per i 6= j (2.37)

Ψ = 0 se rij ≤ a per i 6= j (2.38)

dove la funzione d’onda Ψ = Ψ(~r1, · · · , ~rN) è funzione di tutte le coordinate.L’interazione di tipo sfera rigida è equivalente ad una condizione al contorno

15

Figura 2.1: Funzioni d’onda per un potenziale con lunghezza di scatteringpositiva (a sinistra) o negativa (a destra) Figura tratta da [1] (fig. 10.3 e10.4).

che richieda l’annullamento della funzione d’onda per rij ≤ a per tutti glii 6= j. Nello spazio delle configurazioni 3N -dimensionale l’insieme dei puntiche per cui rij = a rappresenta una ipersuperficie della forma di un “cilindro”(con le coordinate diverse da ~ri e ~rj libere di variare); in tutto ci sono N(N−1)/2 di questi “cilindri” che si intersecano in modo complicato dando luogoad una struttura ad “albero” (Fig. 2.2).

Analogamente a quanto fatto prima, si può sostituire l’effetto di ognicilindro con multipoli lungo gli assi equivale ad introdurre degli pseudopo-tenziali a due corpi descritti nel paragrafo precedente. Tuttavia questi nondanno il corretto andamento di Ψ in prossimità dell’intersezione di due o piùcilindri, che corrispondono ad urti a più corpi di cui non si tiene conto deglipseudopotenziali a due corpi. Per avere il giusto comportamento anche nellezone di intersezione è quindi necessario introdurre degli pseudopotenziali atre o più corpi che devono essere ricavati risolvendo il problema corrispon-dente. La dipendenza di questi pseudopotenziali da a, diametro della sferarigida, può però essere facilmente ricavata da considerazioni dimensionali.

Ad esempio lo pseudopotenziale a tre corpi deve apparire nell’equazionedi Schödinger a tre corpi nella forma

(∇21 +∇22 +∇23 + k2)Ψ=∑

(pseudopotenziali a 2 corpi) + δ(~r1 − ~r2)δ(~r1 − ~r3)KΨ

La quantità K dovrà avere le dimensioni di una lunghezza alla quarta po-tenza. A basse energie, cioè per k → 0, l’unica lunghezza del problema è a equindi K deve essere dell’ordine di a4. Allo stesso modo si possono trovarele dimensioni degli pseudopotenziali di ordine superiore.

16

Figura 2.2: Rappresentazione schematica dello spazio delle configurazione a3N dimensioni. L’interazione tra sfere rigide è equivalente a condizioni alcontorno nulle per la funzione d’onda sulla superficie dei cilindri. Figuratratta da [1] (fig. 10.5).

Queste considerazioni continuano a valere anche se il potenziale tra par-ticelle non è quello della sfera rigida, ma un potenziale generico a raggiod’azione finito che non abbia stati legati. Inoltre gli pseudopotenziali checoinvolgono più di due particelle possono essere ignorati se si è interessati adun’accuratezza all’ordine a (per un potenziale interatomico generico, a2 perle sfere rigide).

In definitiva quindi (ricordando la 2.24 e tenendo conto che la massaridotta del sistema a due corpi è µ = m/2), l’hamiltoniana efficace per ungas interagente di N particelle identiche di massa m si può scrivere come

H = − ~2

2m(∇21 + · · ·+∇2N) +

4πa~2

m

∑i

della delta di Dirac), abbiamo(∂

∂rrψ

)∣∣∣∣r=0

= ψ(0) + r

(∂ψ

∂r

)∣∣∣∣r=0

= ψ(0) (2.40)

e possiamo quindi trascurare il termine di derivata

∂

∂rijrij ∼ 1 (2.41)

Questo è il caso che si ha per esempio quando l’interazione produce unapiccola perturbazione rispetto alle funzione d’onda del sistema libero, che èregolare in r = 0.

Nel seguito noi supporremo che queste condizioni siano soddisfatte e uti-lizzeremo quindi come hamiltoniana di riferimento per il sistema ad N corpila seguente espressione

H = − ~2

2m(∇21 + · · ·+∇2N) +

4πa~2

m

∑i

Capitolo 3

La teoria di Gross-Pitaevskii

In questo capitolo partiremo dai risultati dei capitoli precedenti per introdur-re la teoria che descrive i condensati di Bose-Einstein. In particolare quindiconsideriamo un sistema di N bosoni con interazioni a due corpi, nell’ipotesidi gas diluito ((N/V )a3 � 1), in cui il potenziale di interazione assume laforma

v(~r − ~r ′) = gδ(~r − ~r ′) con g = 4π~2a

m(3.1)

dove a è la lunghezza di scattering e m la massa degli atomi del gas. Nellaformulazione della teoria dei campi, l’hamiltoniana che descrive il nostrosistema può essere scritta cos̀ı:

Ĥ =

∫d~r ψ̂†(~r)

[− ~

2

2m∇2 + V (~r)

]ψ̂(~r) +

g

2

∫d~r ψ̂†(~r)ψ̂†(~r)ψ̂(~r)ψ̂(~r). (3.2)

L’evoluzione temporale dell’operatore di campo è determinata dall’equa-zione di Heisenberg:

i~∂

∂tψ̂(~r, t) =

[ψ̂(~r), Ĥ

](3.3)

che può essere messa in forma esplicita calcolando il commutatore (ricordan-do la (1.5)). Per quanto riguarda il contributo del termine cinetico e delpotenziale, si ottiene [

− ~2

2m∇2 + V (~r)

]ψ̂(~r) (3.4)

mentre per la parte delle interazioni si ha[ψ̂(~r),

g

2

∫d~r ′ ψ̂†

2

(~r ′)ψ̂2(~r ′)

]=

=g

2

∫d~r ′

[ψ̂(~r), ψ̂†

2

(~r ′)]ψ̂2(~r ′) =

= gψ̂†(~r)ψ̂2(~r). (3.5)

19

Utilizzando le (3.4) e (3.5), l’equazione di Heisenberg (3.3) diventa

i~∂

∂tψ̂(~r, t) =

[− ~

2

2m∇2 + V (~r) + gψ̂†(~r, t)ψ̂(~r, t)

]ψ̂(~r, t) (3.6)

3.1 L’equazione di Gross-Pitaevskii

Arrivati a questo punto passeremo dalla teoria quantistica dei campi ad unateoria semiclassica, ottenuta introducendo il valore medio quantistico (valored’aspettazione sullo stato quantistico del sistema), φ = 〈ψ̂〉. In generalepossiamo scrivere

ψ̂ = 〈ψ̂〉+ δψ̂ (3.7)

dove il termine δψ̂ rappresenta le fluttuazioni quantistiche, e per definizionerisulta 〈δψ̂〉 = 0.

Come sappiamo dalla Meccanica Statistica, un sistema di bosoni a bassetemperature presenta il fenomeno della condensazione di Bose-Einstein, ca-ratterizzato da un’occupazione macroscopica dello stato fondamentale. Latransizione verso questa nuove fase è associata alla comparsa di un parame-tro d’ordine diverso da zero, che nella formulazione statistica della Teoriadei Campi è proprio il valore medio φ del campo che rappresenta anche lafunzione d’onda del condensato. Analogamente a quanto visto preceden-temente, questo è un oggetto che non descrive le particelle del condensatosingolarmente, bens̀ı il comportamento collettivo di tutto il sistema. Il ter-mine di fluttuazione δψ̂ è invece associato alle particelle non condensate.Per temperature molto basse il contributo di queste ultime può, in primaapprossimazione, essere trascurato (approssimazione di campo medio).

In definitiva, effettuando la sostituzione ψ̂(~r, t) → φ(~r, t) nell’equazione(3.6) otteniamo l’equazione di Gross-Pitaevskii (GPE)

i~∂

∂tφ(~r, t) =

[− ~

2

2m∇2 + V (~r) + g |φ(~r, t)|2

]φ(~r, t) (3.8)

che è una equazione di Schrödinger non-lineare (dove il termine non-lineareè legato alle interazioni tra le particelle) che descrive l’evoluzione tempora-le della funzione d’onda del condensato φ(~r, t). Quest’ultima è legata alladensità di probabilità spaziale semplicemente da

ρ(~r, t) = |φ(~r, t)|2 (3.9)

che è una quantità misurabile direttamente negli esperimenti.Ricordando la definizione dell’operatore numero

N̂ =

∫d~r ψ̂†(~r)ψ̂(~r) (3.10)

20

nell’approsimazione di campo medio abbiamo (ricordiamo che 〈δψ̂〉 = 0)

〈N̂〉 =∫d~r〈(φ+ δψ̂

)† (φ+ δψ̂

)〉

=

∫d~r(|φ|2 + 〈δψ̂†δψ̂〉

)≈

∫d~rρ(~r, t) = N (3.11)

dove N è il numero di particelle del condensato, consistentemente con quantoabbiamo discusso fino ad ora. Per comodità la funzione d’onda del conden-sato può anche essere normalizzata ad uno, riscalandola opportunamenteφ→

√Nφ; questo modifica anche il termine di interazione, g → Ng. Questa

scelta ci sarà utile in seguito.L’equazione di Gross-Pitaevskii può anche essere ricavata tramite un

principio variazionale

i~∂

∂tφ(~r, t) =

∂E

∂φ∗(3.12)

dal funzionale energia E

E[φ] =

∫d~r

{φ∗(~r, t)

[− ~

2

2m∇2 + V (~r)

]φ(~r, t) +

g

2|φ(~r, t)|4

}(3.13)

=

∫d~r

{~2

2m|∇φ(~r, t)|2 + V (~r) |φ(~r, t)|2 + g

2|φ(~r, t)|4

}(3.14)

dove il primo termine rappresenta l’energia cinetica Ekin, il secondo l’energiapotenziale Epot, ed il terzo l’energia di interazione Eint.

Si definisce anche la quantità Erel = Ekin + Eint che prende il nome direlease energy (energia di rilascio), che corrisponde all’energia disponibiledurante l’espansione in seguito allo spegnimento del potenziale di trappo-la (questo è un ingrediente essenziale negli esperimenti dove l’imaging delcondensato viene fatto dopo un certo tempo di volo, o time-of-flight).

3.2 L’equazione di Gross-Pitaevskii staziona-

ria

Passiamo ora a discutere lo stato fondamentale del sistema. Questo si ottienecercando una soluzione stazionaria della forma:

φ(~r, t) = φ(~r)e−iµt/~ (3.15)

che, sostituita nell’equazione (3.8), da luogo all’equazione di Gross-Pitaevskiiindipendente dal tempo[

− ~2

2m∇2 + V (~r) + g |φ(~r)|2

]φ(~r) = µφ(~r) (3.16)

21

In questa equazione µ è il potenziale chimico, cioè l’energia necessaria per ag-giungere un atomo al condensato, ed è diverso dall’energia totale del sistema.Moltiplicando a sinistra per φ∗ e integrando in d~r si ottiene infatti

Ekin + Epot + 2Eint = µN (3.17)

da cui

µ =Ekin + Epot + 2Eint

N=E

N+EintN

(3.18)

doveE = Ekin + Epot + Eint (3.19)

è l’energia totale definita sopra. Il potenziale chimico µ quindi differiscedall’energia totale per particella E/N per un termine pari all’energia di inte-razione per particella, e questa differenza è dovuta al fatto che l’equazione diGross-Pitaevskii è non lineare. (Nel caso non interagente, g = 0, l’equazionedi Gross-Pitaevskii si riduce all’equazione di Schrödinger lineare).

Osserviamo anche che la dipendenza di E da N può essere resa esplicitariscalando opportunamente il parametro d’ordine come abbiamo accennatoin precedenza, φ→

√Nφ

E[φ] = N

∫d~r

{~2

2m|∇φ|2 + V (~r) |φ|2 +N2 g

2

∫d~r |φ|4

}≡ NE0 +N2I (3.20)

da cui è facile verificare (dalla (3.18) si ha µ = E0 + 2NI) che vale

∂E

∂N= µ (3.21)

consistentemente con l’interpretazione termodinamica di µ come energia ne-cessaria per aggiungere una particella al condensato.

Facciamo ora qualche considerazione su g. Innanzitutto ricordiamo che

g =4π~2am

(3.22)

A seconda del segno della lunghezza di scattering a si ha una interazioneattrattiva (a < 0) o repulsiva (a > 0). In particolare, per a < 0 l’interazionetende a concentrare tutti gli atomi sono in una regione limitata di spazio,sempre più piccola quanto più è forte l’attrazione. Questo è possibile fino adun certo valore critico, che identifica la condizione per cui non esiste più unasoluzione normalizzabile dell’equazione (3.16). Dal punto di vista fisico, inquesta situazione anche le interazioni a tre corpi diventano importanti ed acausa di esse si ha una perdita di atomi del condensato.

Nel seguito della discussione considereremo sempre il caso a > 0.

22

3.2.1 Il potenziale armonico

Uno dei potenziali più tipici utilizzati negli esperimenti con i condensati è unpotenziale armonico del tipo

Vho =1

2m∑i

ω2i r2i (3.23)

ed è inoltre utile definire ωho = 3√ωxωyωz.

Nel caso il cui il sistema ha una simmetria cilindrica (ωx = ωy ≡ ω⊥) ilpotenziale può essere scritto come

Vho =1

2mω2⊥r

2⊥ +

1

2mω2zz

2 (3.24)

mentre nel caso di simmetria sferica si semplifica ulteriormente ponendo ω⊥ =ωz = ωho

Vho =1

2mω2hor

2 . (3.25)

Per semplificare la risoluzione dell’equazione di Gross-Pitaevskii è possi-bile passare a variabili adimensionali; consideriamo per semplicità il caso disimmetria sferica e definiamo

aho =

√~

mωho(3.26)

che rappresenta la lunghezza caratteristica dell’oscillatore armonico. Intro-duciamo quindi la variabile adimensionale r̃ = r/aho e riscriviamo l’equazione(3.16) come segue[

− ~2

2m

1

a2ho∇2r̃ +

1

2mω2hoa

2hor̃

2 +4π~2am

N

a3ho|φ̃|2]φ̃ = µφ̃ (3.27)

dove la φ̃ è tale che

|φ̃|2 = a3ho

N|φ|2 (3.28)

in modo da essere normalizzata ad uno.

Semplificando e dividendo entrambi i termini per ~ωho, equazione diGross-Pitaevskii stazionaria assume la seguente forma adimensionale[

−12∇2r̃ +

1

2r̃2 + 4π

Na

aho|φ̃|2]φ̃ = µ̃φ̃ (3.29)

con µ̃ = µ/~ωho.

23

3.2.2 Il limite Thomas-Fermi

Nel caso in cui il termine di interazione sia predominate su quello cinetico- il cosiddetto limite di Thomas-Fermi - l’equazione di Gross-Pitaevskii sta-zionaria (3.16) può essere ulteriormente semplificata trascurando appunto iltermine cinetico [

V (~r) + g |φ(~r)|2]φ(~r) = µφ(~r) (3.30)

da cui si può ricavare la forma della distribuzione spaziale del condensato

|φ|2 = g−1 (µ− Vho(~r)) =µ

g

(1−

∑i

mω2i2µ

r2i

)(3.31)

che è l’equazione di un paraboloide. Conviene inoltre definire il raggio diThomas-Fermi nella direzione i-esima come

RTFi =

√2µ

mω2i(3.32)

per cui la distribuzione di densità assume la forma

|φ|2 = µg

(1−

∑i

(riRTFi

)2). (3.33)

In Fig. 3.1 è mostrato come il profilo della funzione d’onda si modificapassando dal caso non interagente (a = 0, profilo gaussiano) al limite TFdove la forte interazione repulsiva tra gli atomi provoca un allargamento delcondensato, caratterizzato da un profilo parabolico. In realtà il profilo non èperfettamente parabolico, in quanto la distribuzione mantiene delle code.

Nel limite TF il potenziale chimico si ottiene analiticamente dalla condi-zione di normalizzazione, come mostrato qui sotto.

N =

∫d3r|φ|2 = µ

g

∫d3r

(1−

∑i

(riRTFi

)2)(3.34)

=µ

g

∏i

RTFi

∫d3ξ

(1−

∑i

ξ2i

)(3.35)

=

(2µ

~ωho

)2aho15a

(3.36)

Nel limite TF abbiamo quindi

µ =~ωho

2

(15N

a

aho

)2/5∝ N2/5 (3.37)

24

Figura 3.1: Profilo di densità del condensato nel passare dal caso non inte-ragente (a = 0, linea tratteggiata) al limite Thomas-Fermi (a � 1). Figuratratta da [1] (fig. 9).

da cui, ricordando che µ = ∂E/∂N , si ha

E =

∫µdN =

5

7Nµ (3.38)

ed inoltre (Eint = µN − E)

Eint = Erel =2

7Nµ (3.39)

Healing Length

Un’altra lunghezza caratteristica del condensato è la cosiddetta healing lengthξ. Questa si ottiene considerando un condensato con densità uniforme n, edimponendo che la densità vada a zero in un punto arbitrario (questa è peresempio la situazione che si ha ai bordi di un contenitore). La healing lengthrappresenta la distanza su cui la densità ritorna da zero al valore uniforme n, epuò essere ottenuta considerando l’equazione di Gross-Pitaevskii stazionaria[

− ~2

2m∇2 + g |φ(~r)|2

]φ(~r) = µφ(~r) (3.40)

e comparando il termine di energia cinetica e quello di interazione. Per unsistema uniforme la soluzione di questa equazione è φ =

√n, con potenziale

chimico µ = gn (a cui contribuisce solo l’energia di interazione perché quellacinetica è nulla). Quando introduciamo gli effetti ai bordi, dove il termine di

25

interazione si annulla, l’energia mancante deve essere compensata dal terminecinetico. Possiamo quindi scrivere

− ~2

2m∇2 ≈ ~

2

2mξ2= gn =

4π~2am

n (3.41)

da cui risulta

ξ =1√

8πna. (3.42)

La healing length rappresenta la scala di lunghezza caratteristica associataal campo medio (il termine di interazione).

3.3 Formulazione idrodinamica

In questa sezione vedremo che l’equazione di Gross-Pitaevskii può essere ri-condotta a due equazioni ben note, cioè l’equazione di continuità e l’equazionedi Eulero per un superfluido. Per ricavarle, possiamo riscrivere la funzioned’onda del condensato nella rappresentazione modulo-fase

φ =√neiS (3.43)

dove la densità n e la fase S sono entrambe funzioni di ~r e t, n = n(~r, t) eS = S(~r, t). In particolare la fase è legata alla velocità. Infatti

~j = n~v =~

2im(φ∗∇φ− φ∇φ∗) =

=~

2im

[√ne−iS

(∇√n)eiS +

√ni(∇S)eiS

)− c.c.

]=

~2im

[√n(∇

√n) + ni(∇S)− c.c.

]=

~mn∇S (3.44)

da cui semplicemente:

~v =~m∇S (3.45)

Vediamo ora come si trasforma l’equazione di Gross-Pitaevskii dipendentedal tempo in eq. (3.8), che riscriviamo qui sotto

i~∂

∂tφ =

~2

2m∇2φ+ V φ+ g|φ|2φ; (3.46)

utilizzando la forma di φ in (3.43), il primo membro diventa

i~(∂

∂t

√n

)eiS − ~

√n

(∂

∂tS

)eiS, (3.47)

26

mentre per il secondo membro abbiamo

− ~2

2m∇[(∇√n)eiS + i (∇S)

√neiS

]= − ~

2

2meiS[(∇2√n)

(3.48)

+i (∇S)(∇√n)

+ i(∇2S

)√n+ i (∇S)

(∇√n)− (∇S)2

√n].

Gli altri due termini contengono il potenziale esterno di interazione e ilpotenziale interatomico

V√neiS + gn

√neiS. (3.49)

A questo punto uguagliando le parti reali ed immaginarie dei due membridell’equazione di Gross-Pitaevskii, si ottengono due equazioni che devonoessere soddisfatte contemporaneamente

~∂

∂t

√n = − ~

2

2m

(2(∇√n)∇S +

√n∇2S

)(3.50)

−~√n∂

∂tS = − ~

2

2m

(∇2√n−

√n(∇S)2

)+ V

√n+ gn

√n (3.51)

Dalla (3.50), semplificando ~ e moltiplicando entrambi i membri per√n, si

ottiene

√n∂

∂t

√n+

~2m

(2√n(∇√n)∇S + n∇2S

)= 0

1

2

∂

∂t

(√n√n)

+~

2m

(21

2

(∇√n√n)∇S + n∇2S

)= 0

∂n

∂t+

~m

((∇n)∇S + n∇2S

)= 0

∂n

∂t+

~m∇ (n∇S) = 0

∂n

∂t+∇ (n~v) = 0 (3.52)

che è l’equazione di continuità.Dalla (3.51) invece si ha, moltiplicando per

√n

n~∂

∂tS − ~

2

2m

(√n∇2

√n− n(∇S)2

)+ V n+ gn2 = 0

~∂

∂tS − ~

2

2m

(√n (∇2

√n)

n− (∇S)2

)+ V + gn = 0

e facendone il gradiente

~∂

∂t∇S +∇

[− ~

2

2m

√n (∇2

√n)

n+

~2

2m(∇S)2 + V + gn

]= 0

27

ricordando che ∇S = (m/~)~v si ottiene infine

∂

∂t~v +

1

m∇[− ~

2

2m

(∇2√n)√n

+1

2m~v2 + V + gn

]= 0 (3.53)

che è analoga all’equazione idrodinamica di Eulero di fluido ideale senzaviscosità (un superfluido)

∂

∂t~v = − 1

mn∇p− 1

2∇~v2 − 1

m∇V (3.54)

dove il termine di pressione p contiene il contributo delle interazioni gn2/2 edun contributo puramente quantistico −(~2/2m)(∇2

√n)/

√n, detto pressione

quantistica. Osserviamo che la (3.54) è scritta qui per il caso ∇ × ~v = 0(∇×∇S = 0).

Riassumendo quindi, nella formulazione idrodinamica l’equazione di Gross-Pitaevskii (3.8) è equivalente al sistema di equazioni (3.52) e (3.53)

∂n

∂t+∇ (n~v) = 0 (3.55)

m∂

∂t~v +∇

[− ~

2

2m

(∇2√n)√n

+1

2m~v2 + V + gn

]= 0 (3.56)

Bibliografia

[1] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, S. Stringari, Theory of Bose-Einstein condensation in trapped gases, Rev. Mod. Phys. 71, 463 (1999).

[2] C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases(Cambridge University Press, Cambridge, 2008).

28

Capitolo 4

Equazioni di scaling

In questo capitolo prenderemo in considerazione la classe di soluzioni dell’e-quazione di Gross-Pitaevskii - o equivalentemente delle equazioni idrodina-miche - che presentano un comportamento di scaling. Tali soluzioni corri-spondono ad avere una distribuzione di densità che mantiene la stessa formafunzionale dalle coordinate spaziale, che vengono però riscalate tramite deifattori di scala dipendenti dal tempo. Questo approccio risulta applicabilenel limite non-interagente ed in quello opposto di Thomas-Fermi, nel casodi confinamento armonico con frequenze dipendenti dal tempo, ed è utileper studiare le piccole oscillazioni (ottenute modulando o cambiando istan-taneamente le frequenze di confinamento) o per l’espansione libera (assenzadi confinamento).

4.1 Scaling delle equazioni idrodinamiche

Iniziamo considerando il caso di un condensato confinato in un potenzialearmonico dipendente dal tempo, della forma

V (~r, t) =1

2m

3∑i=1

ω2i (t)r2i ; (4.1)

supporremo anche che all’istante iniziale il condensato abbia un profilo Thomas-Fermi

n0(~r) =µ− V0(~r)

g. (4.2)

In questo caso, come verificheremo in seguito, una buona ipotesi di scalingper la densità corrisponde a richiedere

n(~r, t) =1∏j λj

n0

(riλi(t)

)(4.3)

29

in cui la dipendenza dal tempo è contenuta solo nei coefficienti λi(t) cheriscalano le coordinate (il fattore contenente la produttoria serve a conser-vare la normalizzazione). Inseriamo adesso questa ipotesi nell’equazione dicontinuità (3.52); per il termine di derivata temporale si ha

∂n

∂t=∑j

(− λ̇jλj

1∏i λi

n0 +1∏i λi

∂n0∂(rj/λj)

rjλ2j

(−λ̇j)

)

= − 1∏i λi

∑j

λ̇jλj

(n0 + rj

∂n0∂rj

)(4.4)

mentre per il termine di gradiente

∇(n~v) = ∇n · ~v + n∇ · ~v =∑j

[(∇jn)vj + n(∇jvj)] =

=1∏i λi

∑j

[(∇jn0)vj + n0(∇jvj)] (4.5)

dove il simbolo ∇j indica la derivata rispetto a rj.Globalmente, dall’equazione di continuità si ottiene quindi

1∏i λi

∑j

[(− λ̇jλj

+∇jvj

)n0 +

(− λ̇jλjrj + vj

)∇jn0

]= 0 (4.6)

che è soddisfatta se il campo di velocità soddisfa il seguente scaling

vj(~r, t) =λ̇jλjrj; (4.7)

questa equazione completa l’ipotesi di scaling. Osserviamo che fino ad oranon abbiamo utilizzato né l’ipotesi di potenziale armonico (4.1) né quella delprofilo Thomas-Fermi (4.2). Lo scaling della velocità risulta solo dall’ipotesidi scaling sulla densità e dall’equazione di continuità.

Vediamo ora che cosa si ricava dall’equazione di Eulero (3.54); per iltermine di derivata temporale abbiamo

m∂~v

∂t= mrj

(λ̈jλj−λ̇2jλ2j

)(4.8)

mentre per quello di energia cinetica

∇j1

2m~v2 = m~v∇j~v = m

∑k

vk∇jvk

= m∑k

xkλ̇kλkδkj

λ̇kλk

=

= mrjλ̇2jλ2j

(4.9)

30

Per gli altri due termini occorre utilizzare esplicitamente sia l’ipotesi dipotenziale armonico (4.1)

∇V (~r) = mω2j rj (4.10)

che quella del profilo Thomas-Fermi (4.2)

g∇n = g∏i λi∇jn0

(riλi

)=

=g∏i λi

1

g∇j

(µ− 1

2m∑k

ω2k(0)x2kλ2k

)=

= − mrj∏i λi

ω2j (0)

λ2j. (4.11)

Dalla (3.54) si ha quindi

mrj

(λ̈jλj−λ̇2jλ2j

+λ̇2jλ2j−

ω2j (0)

λ2j∏

i λi+ ω2j (t)

)= 0 (4.12)

da cui si ottiene l’equazione per i coefficienti di scaling λj(t)

λ̈j(t) =ω2j (0)

λj∏

i λi− ω2j (t)λj (4.13)

Nota: soluzioni di scaling possono essere ottenute anche in situazioni piùcomplesse, per un potenziale quadratico del tipo:

V (~r, t) = mω2ho

[1

2

∑i,j

riUij(t)rj −∑i

ui(t)ri + u0(t)

]. (4.14)

In questo caso occorre fare un’ipotesi analoga per lo scaling della densitàn(~r, t) e della fase S(~r, t) (o equivalentemente della velocità v(~r, t)), scrivendo

n(~r, t) ∝ 12−∑i,j

riAij(t)rj +∑i

ai(t)ri + n0(t) (4.15)

S(~r, t) ∝ −12

∑i,j

riBij(t)rj +∑i

bi(t)ri + s0(t). (4.16)

31

4.2 Scaling dell’equazione di Gross-Pitaevskii

Le equazioni di scaling possono anche essere ricavate direttamente dall’equa-zione di Gross-Pitaevskii. Nel caso appena considerato è naturale considerarela seguente espressione per la funzione d’onda

ψ(~r, t) =1√∏j λj(t)

ψ0

(riλi(t)

, t

)e

im

2~∑j

λ̇j(t)

λj(t)r2j − iβ(t)

(4.17)

dove in ψ0 abbiamo lasciato una dipendenza esplicita da t ed abbiamo sup-posto uno scaling delle coordinate che deriva direttamente da quello per ladensità, mentre per il termine di fase lo scaling è quello del campo di ve-locità che abbiamo ricavato prima, dove si è tenuto conto della relazione~v = (~/m)∇S(~r). Questo termine di fase definisce and una trasformazionedi gauge locale (dipende dalle coordinate spaziali) che corrisponde al seguen-te shift dell’operatore di impulso P̂ = −i~∇, a cui viene sottratto l’impulsolocale (che si avrebbe anche per un gas classico [2])

P̂ → P̂ +m~v (4.18)

Infatti

−i~∇ψ0eiS = (−i~∇ψ0)eiS + ψ0(~∇S)eiS

= eiS(−i~∇+ ~∇S)ψ0 = eiS(−i~∇+m~v)ψ0 (4.19)

e quindi, per ψ0 → ψ0eiS

〈P̂ 〉ψ0 → 〈P̂ +m~v〉ψ0 . (4.20)

Osserviamo che nella (4.17) compare anche una fase β(t) dipendente so-lo dal tempo, da fissare in modo opportuno tramite l’equazione di Gross-Pitaevskii. Procediamo quindi al calcolo esplicito, inserendo la (4.17) nella(3.8). Per il termine di derivata rispetto al tempo abbiamo (utilizzandonotazioni sintetiche per brevità)

i~∂tψ =

(i~∂t

1√∏λ

)ψ0e

iS−iβ +1√∏λ

(i~Dtψ0) eiS−iβ

−~ 1√∏λψ0

(Ṡ − β̇

)eiS−iβ

=

[(i~∑j

−12

λ̇j∏

k 6=j λk

(∏λ)3/2

)ψ0 +

1√∏λ

(−i~

∑j

λ̇jλ2jrj∇̃jψ0 + i~∂tψ0

)

−~ 1√∏λψ0

(m

2~∑j

(λ̈jλj−λ̇2jλ2j

)r2j − β̇

)]eiS−iβ (4.21)

32

mentre per il laplaciano

− ~2

2m∇2ψ = − ~

2

2m

1√∏λ∇{(∇ψ0 + ψ0i∇S) eiS−iβ

}− ~

2

2m

1√∏λ

{∇2ψ0 + 2i∇ψ0 · ∇S + ψ0i∇2S − ψ0(∇S)2

}eiS−iβ

− ~2

2m

1√∏λ

{(∑j

1

λ2j∇̃2jψ0

)+ 2i

∑j

1

λj

(∇̃jψ0

) m~λ̇jλjrj

+iψ0m

~∑j

λ̇jλj− ψ0

∑j

(m

~λ̇jλjrj

)2 eiS−iβ=

1√∏λ

{− ~

2

2m

∑j

1

λ2j∇̃2jψ0 − i~

∑j

1

λj

(∇̃jψ0

) λ̇jλjrj

−iψ0~2

∑j

λ̇jλj

+ ψ0m

2

∑j

λ̇2jλ2jr2j

}eiS−iβ (4.22)

ed infine

V (~r) + g |ψ(~r, t)|2 = 12m∑j

ω2j (t)r2j +

g∏λ|ψ0|2 (4.23)

L’equazione di Gross-Pitaevskii può essere quindi riscritta cos̀ı (semplifi-cando un fattore comune exp{i(s− β)}/

√∏λ)

i~∂tψ0 = i~∑j

λ̇jλ2jrj∇̃jψ0

+i~2

∑j

λ̇j∏

k 6=j λk∏λ

ψ0 +m

2

∑j

(λ̈jλj−λ̇2jλ2j

)r2jψ0 − ~β̇ψ0

− ~2

2m

∑j

1

λ2j∇̃2jψ0 − i~

∑j

1

λj

(∇̃jψ0

) λ̇jλjrj − i

~2

∑j

λ̇jλjψ0 +

m

2

∑j

λ̇2jλ2jr2jψ0

+1

2m∑j

ω2j (t)r2jψ0 +

g∏λ|ψ0|2ψ0

= − ~2

2m

∑j

1

λ2j∇̃jψ0 +

1

2m∑j

(ω2j (t) +

λ̈jλj

)r2jψ0 +

g∏λ|ψ0|2ψ0 − ~β̇ψ0

Osserviamo che fino a questo punto abbiamo sfruttato solo la forma delpotenziale in eq. (4.1) ma non abbiamo detto niente sulla condizione inizialeψ0. Utilizzando l’ultima equazione è possibile ricavare i comportamenti discaling sia nel caso non-interagente che nel limite Thomas-Fermi.

33

4.2.1 Scaling Thomas-Fermi

Utilizzando l’equazione di scaling (4.13) già trovata in precedenza, e ponendo

~β̇ =µ∏j λj

(4.24)

si ha

i~∂tψ0 = −~2

2m

∑j

1

λ2j∇̃jψ0 +

1

2m∑j

ω2j (0)

λ2j∏λr2jψ0 +

g∏λ|ψ0|2ψ0 −

µ∏λψ0

= − ~2

2m

∑j

1

λ2j∇̃jψ0 +

1∏λ

[1

2m∑j

ω2j (0)r̃2j + g|ψ0|2 − µ

]ψ0

= − ~2

2m

∑j

1

λ2j∇̃jψ0 +

1∏λ

[V (~̃r, 0) + g|ψ0|2 − µ

]ψ0 (4.25)

per cui, se siamo nel limite Thomas-Fermi ed il condensato si trova inizial-mente nel ground-state (4.2), il termine di derivata spaziale è trascurabile,e il termine tra parentesi quadra si annulla per la (4.2). Questo implicai~∂tψ0 = 0, ovvero che non c’è una dipendenza esplicita dal tempo, che ètutta contenuta nei parametri di scaling. Se invece siamo prossimi ma nonesattamente nel limite TF, l’equazione (4.25) è comunque utile per risolverela dinamica residua che non è già contenuta nello scaling1.

4.2.2 Scaling non-interagente

Per procedere in modo analogo al caso precedente, dobbiamo individuareun’equazione per i parametri di scaling λi(t) che ci permetta di far comparirel’equazione stazionaria per il caso non-interagente; poniamo quindi

~β̇ =∑j

1

2

~ωj(0)λ2j

(4.26)

e

λ̈j(t) =ω2j (0)

λ3j− ω2j (t)λj (4.27)

da cui si ha (g = 0)

i~∂tψ0 = −~2

2m

∑j

1

λ2j∇̃jψ0 +

1

2m∑j

ω2j (0)

λ4jr2jψ0 −

∑j

1

2

~ωj(0)λ2j

ψ0

=∑j

1

λ2j

[− ~

2

2m∇̃j + V (~̃r, 0)−

1

2~ωj(0)

]ψ0 (4.28)

1Nel caso dell’espansione libera, che tratteremo in seguito, questo permette di ridurrenotevolmente le difficoltà numeriche in quanto l’aumento delle dimensioni del sistema - equindi della griglia di calcolo - è compensato dai parametri di scaling.

34

anche questa volta, se il sistema si trova inizialmente nel ground-state (non-interagente), il termine tra parentesi quadre è nullo e tutta la dipendenza daltempo è quindi contenuta nelle equazioni di scaling (4.27). Osserviamo che,a differenza del caso Thomas-Fermi, nel caso non interagente le equazioniper i parametri di scaling λj sono disaccoppiate (non compare il termine diproduttoria al denominatore).

4.3 Soluzione delle equazioni di scaling

Abbiamo visto come dall’equazione di Gross-Pitaevskii possiamo ricavare duediverse equazioni di scaling per l’evoluzioni dei parametri λj(t) nel caso dellimite di Thomas-Fermi e nel caso non interagente

λ̈j(t) =ω2j (0)

λj∏

i λi− ω2j (t)λj limite Thomas-Fermi

λ̈j(t) =ω2j (0)

λ3j− ω2j (t)λj limite non interagente

4.3.1 Espansione libera

Come primo esempio dell’utilizzo delle equazioni di scaling consideriamo l’e-spansione libera di un condensato inizialmente confinato da un potenzialearmonico. Supporremo quindi che all’istante t = 0 le frequenze di trappolasiano poste a zero

ωi(t) =

{ωi(0) 6= 0 t ≤ 00 t >0 .

Supporremo inoltre che il condensato si trovi nel ground-state, con condizioniiniziali che quindi sono λj(0) = 1 e λ̇j(0) = 0.

Caso non interagente. Con queste ipotesi l’evoluzione dei parametridi scaling è data dall’equazione

λ̈j(t) =ω2j (0)

λ3j(4.29)

la cui soluzione è

λj(t) =√

1 + ω2j (0)t2 (4.30)

come si può facilmente verificare.

35

Dimostrazione.

λ̇j(t) =1

2

2ω2j (0)t√1 + ω2j (0)t

2(4.31)

λ̈j(t) =ω2j (0)√

1 + ω2j (0)t2− 1

2

ω2j (0)t2ω2j (0)t(√

1 + ω2j (0)t2)3

=1

λ3j

[ω2j (0)

(1 + ω2j (0)t

2)− ω4j (0)t2

]=ω2j (0)

λ3j(4.32)

Questo vuol dire che la larghezza del condensato, che è inizialmente fis-sata2 dalla lunghezza dell’oscillatore in ciascuna direzione Rj ∝

√~/mωj,

evolve nel tempo come

Ri(t) = Ri(0)√

1 + ω2i (0)t2

∝√

~mωi

√1 + ω2i (0)t

2 (4.33)

Consideriamo in particolare il caso di un potenziale a simmetria cilindrica

V =1

2mω2⊥r

2⊥ +

1

2mω2zz

2 ; (4.34)

in questo caso si usa definire l’aspect ratio come rapporto tra i raggi quadraticimedi radiale e assiale

R⊥(t)

Rz(t)=

√ωzω⊥

√1 + ω2⊥t

2

1 + ω2zt2

=√λ

√1 + ω2⊥t

2

1 + ω2zt2

(4.35)

con λ = ωz/ω⊥. È facile quindi verificare che mentre l’aspect ratio iniziale è√λ, nel limite t→∞ (t� 1/ωz, 1/ω⊥) si ha

R⊥(t)

Rz(t)−−−→t→∞

√ωrωz

=1√λ

(4.36)

In sostanza quello che accade è che nel limite di t→∞ l’aspect ratio divental’inverso di quello iniziale (vedi linea tratteggiata in Fig. 4.1).

2Nel caso non interagente la “larghezza” del condensato - che si trova nel ground-statedell’oscillatore armonico - è data dal raggio quadratico medio Ri =

√〈r2i 〉. Poiché l’unica

scala di lunghezze è quella fissata dalla lunghezza dell’oscillatore ahoj =√

~/mωj , saràRi ∝

√~/mωj , con una costante di proporzionalità che è un semplice fattore numerico.

36

Figura 4.1: Evoluzione dell’aspect ratio ottenuta dalla soluzione delle equa-zioni di scaling per l’espansione libera, nel caso non interagente (linea tratteg-giata) e in quello Thomas-fermi (linea continua), per due diverse situazionisperimentali (rappresentate dai punti, ben descritti dal regime TF). Figuratratta da [1] (fig. 18).

Limite Thomas-Fermi. In questo caso la presenza di un termine diaccoppiamento tra i diversi parametri di scaling λj rende più difficile unarisoluzione generale del problema. Considereremo quindi il caso particolaredel potenziale a simmetria cilindrica che abbiamo appena visto per il casonon interagente. Con questa ipotesi il sistema in Eq. (4.13) si riduce a duesole equazioni

λ̈⊥ =ω2⊥(0)

λ3⊥λz(4.37)

λ̈z =ω2z(0)

λ2zλ2⊥

(4.38)

Supporremo inoltre di essere in geometria fortemente allungata nella direzio-ne assiale, ωz(0)/ω⊥(0) ≡ �� 1. Passando quindi a variabili adimensionali,ponendo τ = ω⊥(0)t, le equazioni precedenti si riscrivono cos̀ı

d2λ⊥dτ 2

=1

λ3⊥λz(4.39)

d2λzdτ 2

=�2

λ2zλ2⊥

(4.40)

con condizioni iniziali analoghe al caso precedente, λj(0) = 1, λ̇j(0) = 0.

37

Cerchiamo ora una soluzione che abbia la forma di uno sviluppo in serie di �

λ⊥(t) =∞∑n=0

�nλ(n)⊥ (t) (4.41)

λz(t) =∞∑n=0

�nλ(n)z (t) (4.42)

Riscriviamo il sistema di equazioni differenziali in forma compatta{λ3⊥λzλ̈⊥ = 1

λ2⊥λ2zλ̈z = �

2

sviluppandolo all’ordine più basso in �. Scriviamo quindi la prima equazioneall’ordine zero e la seconda fino al secondo ordine in �

λ(0)⊥

3λ(0)z λ̈

(0)⊥ = 1 (4.43)(

λ(0)⊥ + �λ

(1)⊥ + �

2λ(2)⊥

)2·

·(λ(0)z + �λ

(1)z + �

2λ(2)z)2(

λ̈(0)z + �λ̈(1)z + �

2λ̈(2)z

)2= �2 (4.44)

Consideriamo la (4.44); all’ordine zero questa è semplicemente

λ(0)⊥

2λ(0)z

2λ̈(0)z = 0 (4.45)

poichè λ(0)⊥ e λ

(0)z sono diverse da zero, deve essere

λ̈(0)z = 0 (4.46)

che implica, ricordando le condizioni iniziali, che λ̇(0)z (t) = cost = λ̇

(0)z (0) = 0.

A sua volta questo implica che λ(0)z (t) = cost = λ

(0)z (0) = 1; inoltre la (4.43)

diventa

λ̈(0)⊥ =

1

λ(0)⊥

3 (4.47)

che è analoga all’equazione del caso non interagente e quindi ha come solu-zione

λ(0)⊥ =

√1 + τ 2 =

√1 + ω2⊥(0)t

2 (4.48)

Riscriviamo la (4.44) tenendo conto che λ̈(0)z ≡ 0 e λ(0)z ≡ 1(

λ(0)⊥ + �λ

(1)⊥ + �

2λ(2)⊥

)2(1 + �λ(1)z + �

2λ(2)z)2(

�λ̈(1)z + �2λ̈(2)z

)2= �2 (4.49)

e consideriamo il primo ordine in �

λ(0)⊥

2λ̈(1)z = 0 (4.50)

38

cioè, come per l’equazione all’ordine zero

λ̈(1)z = 0 (4.51)

da cui λ̇(1)z (t) = cost = λ̇

(1)z (0) = 0; di conseguenza anche il termine di ordine

� di λz è nullo.Consideriamo infine il termine di ordine �2 della (4.44). Dato che anche

λ̈(1)z = 0, questa diventa

λ(0)⊥

2�2λ̈(2)z = �

2 (4.52)

da cui

λ̈(2)z =1

λ(0)⊥

2 =1

1 + τ 2(4.53)

la cui soluzione èλ(2)z (τ) = τ arctan τ − ln

√1 + τ 2 (4.54)

come è facile verificare.

Dimostrazione.

dλzdτ

= arctan τ +τ

1 + τ 2− 1√

1 + τ 21

2

1√1 + τ 2

2τ = arctan τ (4.55)

d2λzdτ 2

=1

1 + τ 2(4.56)

In definitiva quindi la soluzione delle equazioni di scaling è data daλ⊥(τ) =√

1 + τ 2 (4.57)

λz(τ) = 1 + �2(τ arctan τ − ln

√1 + τ 2

)(4.58)

in cui τ = ω⊥(0)t.Ricordiamo che nel limite di Thomas-Fermi il condensato ha un profilo

di densità parabolico caratterizzato da un raggio di Thomas-Fermi Rj, cheevolverà quindi nel seguente modo

R⊥(τ) = R⊥(0)λ⊥(τ) (4.59)

Rz(τ) = Rz(0)λz(τ) (4.60)

con Rj(0) =√

2µ/mω2j . Per quanto riguarda l’aspect ratio abbiamo dunque

R⊥(τ)

Rz(τ)=

(ωzω⊥

) √1 + τ 2

1 + �2(τ arctan τ − ln

√1 + τ 2

) == �

√1 + τ 2

1 + �2(τ arctan τ − ln

√1 + τ 2

) (4.61)39

e nel limite di τ →∞ si ottiene semplicemente

R⊥(τ)

Rz(τ)−−−→τ→∞

�τ

�2τπ

2

=2

π�(4.62)

Come si vede in Fig. 4.1 questo comportamento è ben diverso da quellodel caso non interagente, e ben descrive gli esperimenti.

4.3.2 Piccole oscillazioni

In questo paragrafo considereremo le piccole oscillazioni di un condensato in-teragente confinato in una trappola armonica a simmetria cilindrica, nel casoin cui venga portato leggermente fuori dalla posizione di equilibrio. Vedremoche il sistema si mette ad oscillare intorno alla posizione di equilibrio confrequenze che sono determinate dalle frequenze della trappola ottica di confi-namento. Poiché questo è un moto collettivo di tutti gli atomi del condensato(che si muovono in fase), sono dette oscillazioni collettive.

Riprendiamo quindi le equazioni di scaling nel limite di Thomas-Fermiλ̈⊥ =

ω2⊥(0)

λ3⊥λz− ω2⊥(t)λ⊥ (4.63)

λ̈z =ω2z(0)

λ2zλ2⊥− ω2z(t)λ⊥ (4.64)

dove supporremo che le frequenze di trappola siano costanti nel tempo,ωj(t) = ωj(0) ≡ ωj. Analogamente al caso dell’espansione libera visto inprecedenza, queste equazioni possono essere riscritte in forma adimensionaledefinendo � = ωz/ω⊥ e τ = ω⊥t, da cui si ha

d2λ⊥dτ 2

=1

λ3⊥λz− λ⊥ (4.65)

d2λzdτ 2

=�2

λ2zλ2⊥− �2λz (4.66)

Nel regime di piccole oscillazioni possiamo scrivere

λj(τ) = 1 + δλj(τ) con δλj � 1 (4.67)

da cui si ottiened2δλ⊥dτ 2

=1

(1 + δλ⊥)3(1 + δλz)

− (1 + λ⊥) (4.68)

d2δλzdτ 2

=�2

(1 + δλz)2(1 + δλ⊥)

2 − �2(1 + δλz). (4.69)

40

Limitandoci al primo ordine in δλ il precedente sistema di equazioni diventad2δλ⊥dτ 2

' −4δλ⊥ − δλz (4.70)

d2δλzdτ 2

' −�2(2δλ⊥ + 3δλz) (4.71)

e può essere riscritto in forma matriciale ponendo

δλ =

(δλ⊥δλz

)(4.72)

da cui si ottiened2

dτ 2δλ = −Mδλ (4.73)

dove M è la seguente matrice 2× 2

M =

(4 1

2�2 3�2

). (4.74)

L’equazione (4.73) ha soluzioni del tipo

δλ = δλ0e−iωτ (4.75)

con ω che è soluzione dell’equazione agli autovalori

det(M − ω2I

)= 0 (4.76)

ovverosia

det

(4− ω2 1

2�2 3�2 − ω2)

= 0 (4.77)

da cui si ottiene

ω2 =4 + 3�2 ±

√9�4 − 16�2 + 162

. (4.78)

All’ordine più basso in �, le due soluzioni (corrispondenti al ±) sono

ω2+ =4 + 4

2= 4 (4.79)

e

ω2− =4 + 3�2 − 4

√1− �2

2=

3�2 + 2�2

2=

5

2�2. (4.80)

A questo punto possiamo far ricomparire esplicitamente la scala di fre-quenze ω⊥, da cui si ha che le frequenze dei due modi di oscillazione sono

ω+ = 2ω⊥ (4.81)

ω− =

√5

2

ωzω⊥

ω⊥ =

√5

2ωz (4.82)

41

Calcolando anche le ampiezze dei due modi, troveremmo che il primo modo- di frequenza 2ω⊥ - è un’oscillazione che avviene principalmente lungo ladirezione radiale (detto transverse breathing), mentre il secondo - di frequenza(5/2)ωz - avviene principalmente lungo la direzioni assiale (detto quadrupolo).Approfondiremo meglio l’analisi dei moti collettivi nel prossimo capitolo.

Per concludere ricordiamo che questi risultati appena discussi sono sta-ti ottenuti con uno sviluppo al primo ordine in δλ. Come abbiamo vistoin questa approssimazione le equazioni sono lineari e i due modi risultanodisaccoppiati; considerando invece anche gli ordini successivi la situazionesi complica perchè si esce dal regime lineare e compaiono anche termini diaccoppiamento tra modi assiali e radiali.

4.4 Scaling di un gas termico non interagente

Comportamenti di scaling possono essere ottenuti anche in presenza di un gastermico3 non interagente in un potenziale armonico dipendente dal tempo [3].In questo caso occorre prendere in considerazione l’equazione di Heisenberg

i~∂

∂tψ̂ =

[ψ̂, Ĥ

]=∑i

[− ~

2

2m∇2i +

1

2mω2i (t)r

2i

]ψ̂ (4.83)

e, come nel caso dell’equazione di Gross-Pitaevskii, possiamo supporre unoscaling dell’operatore di campo del tipo

ψ̂(~r, t) =1√∏j λj(t)

Φ̂

(riλi(t)

, t

)e

im

2~∑j

λ̇j(t)

λj(t)r2j

(4.84)

con la differenza che non abbiamo incluso esplicitamente il fattore di faseβ(t). Richiedendo quindi che valgano le usuali equazioni di scaling per iparametri λi per il caso non interagente

λ̈i(t) =ω2i (0)

λ3i− ω2i (t)λi (4.85)

si ha (x̃i = ri/λi)

i~∂

∂tΦ̂ =

∑i

1

λ2i

[− ~

2

2m∇2x̃i +

1

2mω2i (0)x̃

2i

]Φ̂. (4.86)

3Cioè non a T = 0 come abbiamo considerato formalmente fino ad ora per uncondensato. La trattazione vale sia per bosoni che per fermioni.

42

Dal momento che il problema è separabile, possiamo fattorizzare l’operatoredi campo come

Φ̂ =∏j

φ̂(x̃j, t) (4.87)

e definire inoltre una nuova variabile temporale τ per ciascuna direzionespaziale tramite la relazione

τ =

∫ t0

dt′1

λ2i (t′)

(4.88)

per cui l’equazione di evoluzione dell’operatore di campo (4.86) può esserescomposta nelle tre componenti

i~∂

∂τφ̂(x̃i, τ) =

[− ~

2

2m∇2x̃i +

1

2mω2i (0)x̃

2i

]φ̂(x̃i, τ). (4.89)

Quindi, con queste trasformazioni, il problema è stato ridotto al caso banaledi evoluzione in un potenziale armonico indipendente dal tempo. Supponen-do di partire da una situazione iniziale di equilibrio, per la distribuzione didensità del gas abbiamo

n(~r, t) = 〈ψ̂†ψ̂〉 = 1∏j λj

〈Φ̂†Φ̂〉 = 1∏j λj

n0

(riλi

)(4.90)

dove n0(~r) è la distribuzione all’istante iniziale t = 0, che non dipendeesplicitamente dal tempo essendo per ipotesi di equilibrio.

4.4.1 Distribuzione di equilibrio

Nell’approssimazione semiclassica, la distribuzione di equilibrio può esserescritta cos̀ı (β = 1/kBT )

n0(~r, ~p) =1

eβ(ε(~r,~p)−µ) ± 1(4.91)

dove il segno più è per la statistica di Fermi-Dirac mentre il segno meno èper la statistica di Bose-Einstein; l’energia ε(~r, ~p) è data da

ε(~r, ~p) =p2

2m+ Vho(~r). (4.92)

La densità spaziale si ottiene quindi semplicemente integrando sugli impulsi

n0(~r) =

∫d3p

(2π~)31

eβ(ε(~r,~p)−µ) ± 1(4.93)

43

da cui

n0(~r) =

1

λ3Tg 3

2

(eβ(µ−Vho(~r))

)per bosoni (4.94)

1

λ3Tf 3

2

(eβ(µ−Vho(~r))

)per fermioni (4.95)

dove λT è la lunghezza d’onda termica

λT =

√2π~2mkBT

e le funzioni g 32

e f 32

sono definite nel seguente modo

gα(z) =∞∑n=1

zn

nα= Liα(z) (4.96)

fα(z) = −Liα(−z) =1

Γ(α)

∫dy

yα−1

(z−1ey + 1), (4.97)

dove Γ(α) è la funzione Gamma di Eulero. Notiamo che l’espressione per ibosoni è vale per temperature maggiori della temperatura di condensazione,avendo omesso il termine relativo alle particelle nel condensato.

A questo punto possiamo calcolare la larghezza quadratica media delladistribuzione spaziale di densità

〈r2i 〉 =∫d3r r2i n0(~r) =

∫d3r

1

λ3T

1

Γ(3/2)

∫dy

y1/2r2i(e−β(µ−Vho(~r))ey ± 1)

=

=1∏

j(√mωj)

1

mω2i

∫d3ξ

λ3T

1

Γ(3/2)

∫dy

y1/2ξ2i(e−β(µ−

Pi

12ξ2i )ey ± 1

) (4.98)dove abbiamo definito ξi =

√mωiri. Abbiamo quindi

〈r2i 〉 ∝1

ω2i∏

j(ωj)(4.99)

dove il fattore di proporzionalità è un fattore numerico, identico per ognidirezione i (dipende anche dalla massa). Inoltre la larghezza quadraticamedia dipende anche inversamente dalla produttoria delle frequenze in tuttele direzioni.

4.4.2 Espansione libera

Consideriamo ora il caso dell’espansione libera. L’evoluzione nel tempo è de-terminata solo dai parametri di scaling, e come nel caso del condensato abbia-mo ri(t) = ri(0)

√1 + ω2i t

2. Se consideriamo quindi l’evoluzione dell’aspect

44

ratio si ottiene √〈r2i (t)〉√〈r2j (t)〉

=ωjωi

√1 + ω2i t

2√1 + ω2j t

2(4.100)

ed è facile vedere che nel limite di t → ∞ l’aspect ratio tende a 1. Nelcaso di un gas termico quindi la distribuzione diventa asintoticamente iso-tropa, a differenza di quello che abbiamo visto per un condensato in cui ladistribuzione si mantiene anisotropa anche asintoticamente.

Approccio semiclassico

Allo stesso risultato si può giungere anche considerando un’espansione bali-stica, cioè un’espansione del sistema in assenza di un potenziale di confina-mento. In questo caso la densità è data dalla relazione

n(~r, ~p, t) = n0

(~r − ~p

mt, ~p

)(4.101)

Consideriamo quindi il termine di energia ε(~r, ~p)

ε(~r, ~p) =p2

2m+∑i

1

2mω2i r

2i →

p2

2m+∑i

1

2mω2i

(ri −

pimt)2

(4.102)

che può essere sviluppato cercando di far comparire p in un unico terminequadratico

p2

2m+∑i

1

2mω2i

(~r − ~p

mt

)2=

=p2

2m+

1

2m∑i

ω2i r2i −

∑i

ω2i ripit+1

2m∑i

ω2ip2i t

2

m2=

=1

2m

∑i

p2i(1 + ω2i t

2)−∑i

ω2i ripit+1

2m∑i

ω2i r2i =

=1

2m

∑i

p2iλ2i −

2m

2m

∑i

ω2i ripit

λiλi +

1

2m∑i

ω2i r2i =

=1

2m

∑i

(piλi −

mω2i rit

λi

)2− 1

2m

∑i

(mω2i rit

λi

)2+

1

2m∑i

ω2i r2i =

=1

2m

∑i

q2i +1

2m∑i

ω2i r2i

(1− ω

2i t

2

λ2i

)≡ ε̄(~r, ~q) (4.103)

45

dove abbiamo definito λi(t) =√

1 + ω2i t2 e qi = piλi −mω2i rit/λi. Semplifi-

cando ulteriormente abbiamo

ε̄(~r, ~q) =1

2m

∑i

q2i +1

2m∑i

ω2i r2i

(1− ω

2i t

2

λ2i

)(4.104)

=1

2m

∑i

q2i +1

2m∑i

ω2ir2iλ2i

(λ2i − ω2i t2

)=

=1

2m

∑i

q2i +1

2m∑i

ω2ir2iλ2i

(4.105)

Possiamo quindi scrivere

n(~r, t) =

∫d3p

(2π~)3n(~r, ~p, t) =

1∏i λi

∫d3q

(2π~)31

eβ(ε̄(~r,~q)−µ) ± 1(4.106)

da cui è facile ottenere

n(~r, t) =1∏j λj

n0

(riλi

). (4.107)

4.5 Distribuzione dei momenti

In quest’ultima sezione vedremo come è possibile ricavare delle informazionisulla distribuzione nello spazio dei momenti di un condensato dalla distribu-zione spaziale dopo un’espansione libera. Nel caso non interagente4, l’evo-luzione libera della funzione d’onda ψ(~r, t) del condensato è descritta dallaseguente equazione di Schrödinger

i~∂

∂tψ = − ~

2

2m∇2ψ (4.108)

la cui soluzione può essere scritta cos̀ı

ψ(~r, t) =

∫d3p ψ̃0(~p)e

i ~p·~r~ e−ip2

2mt~ (4.109)

Consideriamo l’esponente e cerchiamo anche in questo caso di far apparireun termine quadratico; abbiamo

~p · ~x~

− p2

2m

t

~= − t

2m~

(p2 − 2m

t~p · ~x

)=

= − t2m~

(~p− m~x

t

)2+

t

2m~m2x2

t2=

= − t2m~

(~p− m~x

t

)2+

1

2mx2

~t2(4.110)

4Sperimentalmente è possibile annullare il termine nonlineare nell’equazione di Gross-Pitaevskii - che è dovuto alle interazioni a due corpi - tramite dei campi magnetici,sfruttando le cosiddette risonanze di Feshbach che modificano le proprietà collisionali.

46

nel limite t→∞ il secondo termine diventa trascurabile e possiamo scrivere

ψ(~x, t) '∫d3p ψ̃0(~p)e

−i t2m~(~p−

m~xt )

2

=

∫d3q ψ̃0(~q +

m~x

t)e−i

t2m~ q

2

(4.111)

inoltre, poiché nel limite t → ∞ l’esponenziale è un termine rapidamenteoscillante, l’unico termine che conta ai fini dell’integrale è quello che annullal’esponente. Infatti si può dimostrare che vale

limt→+∞

∫dxf(x)e−itx

2

=

√π

teiπ/4f(0) + o(1/t) (4.112)

e di conseguenza, nel limite t→∞ si ha

|ψ(~x, t)|2 ' mht

∣∣∣∣ψ̃0(m~xt)∣∣∣∣2 (4.113)

cioè la distribuzione spaziale nel punto ~x asintoticamente diventa proporzio-nale alla distribuzione iniziale degli impulsi, con ~p = m~x/t. Questa relazioneviene impiegata sperimentalmente - dove il regime asintotico t → ∞ vieneraggiunto dopo qualche decina di millisecondi di espansione balistica - perricostruire la distribuzione degli impulsi. Osserviamo che per ottenere questorisultato non abbiamo fatto ipotesi sul potenziale iniziale.

Bibliografia

[1] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, S. Stringari, Theory of Bose-Einstein condensation in trapped gases, Rev. Mod. Phys. 71, 463 (1999).

[2] Y. Castin and R. Dum, Bose-Einstein condensates in Time DependentTraps, Phys. Rev. Lett. 77, 5315 (1996).

[3] G. M. Bruun and C. W. Clark, Ideal gases in time-dependent traps, Phys.Rev. A 61, 061601 (2000).

47

48

Capitolo 5

Eccitazioni collettive

In questo capitolo discuteremo l’approssimazione di Bogoliubov per descrive-re le eccitazioni di un sistema bosonico, in regime lineare. Questa permettedi calcolare lo spettro delle eccitazioni (fluttuazioni termiche e/o quantisti-che) oltre l’approssimazione di campo medio della teoria di Gross-Pitaevskii.In un sistema uniforme tali eccitazioni risultano essere fononi (onde con re-lazione di dispersione lineare), la cui presenza è strettamente collegata alleproprietà superfluide del sistema. Considereremo infine eccitazioni tipo vor-tici quantizzati (anch’esse manifestazione della superfluidità) per un sistemaconfinato in potenziali armonici (soluzioni dell’equazione di GP non lineare).

5.1 Eccitazioni di Bogoliubov

Consideriamo un gas di bosoni descritto dall’hamiltoniana

Ĥ =

∫d~r ψ̂†(~r)

[− ~

2

2m∇2 + V (~r)

]ψ̂(~r) +

g

2

∫d~r ψ̂†(~r)ψ̂†(~r)ψ̂(~r)ψ̂(~r) (5.1)

dove gli operatori di campo ψ̂ e ψ̂† soddisfano le regole di commutazione[ψ̂(~r), ψ̂†(~r ′)

]= δ(~r − ~r ′)[

ψ̂(~r), ψ̂(~r ′)]

= 0[ψ̂†(~r), ψ̂†(~r ′)

]= 0,

e supponiamo di essere a temperature al di sotto della temperatura critica,dove l’operatore di campo ψ̂ può essere separato nella parte relativa al con-densato ϕ(~r), e in quella delle fluttuazioni (termiche e/o quantistiche, noncondensate) φ̂(~r)

ψ̂(~r) = ϕ(~r) + φ̂(~r) (5.2)

49

dove le regole di commutazione per φ̂(~r) sono le stesse del campo ψ̂(~r), e laϕ soddisfa l’equazione di Gross-Pitaevskii (3.8). Con questa separazione lamatrice densità del sistema diventa

ρ(~r, ~r ′) ≡ 〈ψ̂†(~r)ψ̂(~r ′)〉 = ϕ∗(~r)ϕ(~r) + 〈φ̂†(~r)φ̂(~r ′)〉. (5.3)

A questo punto è utile passare ad una descrizione tipo Meccanica Stati-stica, definendo un’hamiltoniana grancanonica

K̂ = Ĥ − µN̂ = Ĥ − µ∫d~r ψ̂†(~r)ψ̂(~r) (5.4)

dove il potenziale chimico µ, accoppiato all’operatore numero N̂ , ha il ruo-lo di moltiplicatore di Lagrange associato alla conservazione del numero diparticelle (in media).

Utilizzando le (5.1) e (5.2), e tenendo solamente i termini fino al secondoordine in φ, otteniamo la cosiddetta hamiltoniana di Bogoliubov

K̂B =

∫d~r ϕ∗

[− ~

2

2m∇2 + V (~r)− µ+ g

2|ϕ|2

]ϕ

+

∫d~r φ̂†

[− ~

2

2m∇2 + V (~r)− µ+ 2g |ϕ|2

]φ̂

+g

2

∫d~r(ϕ∗2φ̂φ̂+ φ̂†φ̂†ϕ2

)(5.5)

dove abbiamo utilizzato anche il fatto che la funzione d’onda del conden-sato soddisfa l’equazione di Gross-Pitaevskii (3.16). Questo operatore puòessere messo in forma diagonale tramite una trasformazione di Bogoliubov,scrivendo il termine di condensato ϕ come

ϕ =√n0e

iS (5.6)

ed introducendo la seguente trasformazione per i termini φ̂ e φ̂† relativi allefluttuazioni

φ̂(~r) = eiS∑j

[uj(~r)α̂j − v∗j (~r)α̂

†j

](5.7)

φ̂†(~r) = e−iS∑j

[u∗j(~r)α̂

†j − vj(~r)α̂j

](5.8)

dove la somma è estesa solo alla parte non condensata.Vedremo tra poco che gli operatori α̂j e α̂

†j possono essere interpretati co-

me operatori di distruzione e creazioni di quasiparticelle, mentre uj e vj sonodelle funzioni su campo complesso che rappresentano le relative autofunzioni

50

delle eccitazioni al sistema. Gli operatori α̂j e α̂†j devono quindi obbedire alle

regole di commutazione tipiche di operatori di creazione e distruzione[α̂j, α̂

†k

]= δjk (5.9)

[α̂j, α̂k] =[α̂†j, α̂

†k

]= 0. (5.10)

Componendo queste regole di commutazione con le regole di commuta-zione tra φ̂ e φ̂† otteniamo le relazioni di completezza∑

j

[uj(~r)u

∗j(~r

′)− v∗j (~r)vj(~r ′)]

= δ(~r, ~r ′) (5.11)∑j

[uj(~r)v

∗j (~r

′)− v∗j (~r)uj(~r ′)]

= 0 (5.12)∑j

[u∗j(~r)vj(~r

′)− vj(~r)u∗j(~r ′)]

= 0 (5.13)

Sostituendo la (5.7) e la (5.8) nell’hamiltoniana di Bogoliubov (5.5) si ottiene:

K̂B =

∫d~r ϕ∗

[− ~

2

2m∇2 + V (~r)− µ+ g

2|ϕ|2

]ϕ

+∑jk

∫d~r α̂jα̂

†k

[vjLv∗k −

1

2g|ϕ|2 (ujv∗k + u∗kvj)

]+∑jk

∫d~r α̂†jα̂k

[u∗jLuk −

1

2g|ϕ|2

(v∗juk + vku

∗j

)]−∑jk

∫d~r α̂jα̂k

[vjLuk −

1

2g|ϕ|2 (ujuk + vkvj)

]−∑jk

∫d~r α̂†jα̂

†k

[u∗jLv∗k −

1

2g|ϕ|2

(u∗ju

∗k + v

∗kv

∗j

)](5.14)

in cui L è l’operatore hermitiano

L = − ~2

2m(∇+ i∇S)2 − µ− 2g |ϕ|2 (5.15)

e le funzioni u e v obbediscono alle seguenti equazioni{Luj − g |ϕ|2 vj = Ejuj (5.16)Lvj − g |ϕ|2 uj = −Ejvj (5.17)

Si ha inoltre che ∫d~r[u∗j(~r)uk(~r)− v∗j (~r)vk(~r)

]= δij (5.18)∫

d~r [uj(~r)vk(~r)− uk(~r)vj(~r)] = 0 (5.19)

51

Sfruttando le varie proprietà precedenti l’hamiltoniana di Bogoliubov puòessere scritta nella forma diagonale

K̂B =

∫d~r ϕ∗

[− ~

2

2m∇2 + V (~r)− µ+ g

2|ϕ|2

]ϕ

−∑j

∫d~rEj|vj(~r)|2 +

∑j

Ejα̂†jα̂j (5.20)

dove compare esplicitamente il funzionale energia di Gross-Pitaevskii (nellaprima riga) e un termine di hamiltoniana libera (l’ultimo termine) per lequasiparticelle associate agli operatori α̂j, che rappresentano appunto le ec-citazioni termiche/quantistiche (rispettivamente a T 6= 0 e T = 0). Avendoincluso nello sviluppo solo termini fino a quello quadratico, tali eccitazioni ri-sultano non accoppiate tra loro, e si parla quindi di regime lineare (non ci sonotermini non-lineari che le accoppiano). Il termine −

∑j

∫d~rEj |vj(~r)|2, lega-

to alle fluttuazioni quantistiche, prende invece il nome di quantum depletionterm, ed è in prima approssimazione trascurabile.

Per descrivere le proprietà termodinamiche del sistema condensato + ec-citazioni, possiamo quindi utilizzare l’hamiltoniana (5.20), ricordando che inmeccanica statistica il valore di aspettazione di un operatore Ω̂ è dato da

〈Ω̂〉 = TrΩ̂e−βK̂B

Tre−βK̂B(5.21)

da cui ad esempio si ha

〈α̂†jα̂k〉 = 〈α̂kα̂†j〉 − δjk = · · · = δjk

1

eβEj − 1(5.22)

ovverosia che le quasiparticelle seguono la statistica di Bose-Einstein.Abbiamo inoltre

n(~r) = ρ(~r, ~r) = n0(~r) +∑j

|vj(~r)|2 +∑j

1

eβEj − 1(|uj(~r)|2 + |vj(~r)|2

).

(5.23)

5.1.1 Eccitazioni di un condensato uniforme

Calcoliamo ora le eccitazioni di Bogoliubov nel caso di un condensato uni-forme, in assenza di potenziali esterni, V (~r) = 0. In questo caso la soluzionedell’equazione di Gross-Pitaevskii è di tipo onda piana

ϕ =√n0e

i ~p·~r~ =√n0e

i~q·~r (5.24)

e il potenziale chimico è dato da

µ =1

2mv2 + n0g (5.25)

52

Supponiamo che anche le eccitazioni collettive u e v siano del tipo onda piana,cioè

u~k ∼ A~kei~k·~r (5.26)

v~k ∼ B~kei~k·~r (5.27)

Nel caso omogeneo l’operatore L è dato da

L = − ~2

2m(∇+ i~q)2 + 2gn0 −

1

2mv2 − gn0

= − ~2

2m(∇+ i~q)2 + gn0 −

1

2mv2 (5.28)

da cui

Lu = A~kei~k·~r[

~2

2m

(~k + ~q

)2+ gn0 −

1

2mv2

]=

= A~kei~k·~r

[~2k2

2m+

~2~k · ~qm

+~2q2

2m+ gn0 −

1

2mv2

]= A~ke

i~k·~r[ε~k + ~~k · ~v + gn0

](5.29)

dove ~v ≡ ~p/m. In modo analogo possiamo calcolare L∗v, e di conseguenza ilsistema di equazioni {

Lu~k − g |ϕ|2 v~k = E~ku~k

L∗v~k − g |ϕ|2 u~k = −E~kv~k

può essere scritto in forma matriciale

M

(A~kB~k

)= E~k

(A~kB~k

)(5.30)

avendo definito

M =

(ε~k + ~~k · ~v + gn0 −gn0

gn0 −(ε~k − ~~k · ~v + gn0

) ) . (5.31)Gli autovalori si ottengono richiedendo det(M − E~kI) = 0:

−(ε~k + ~~k · ~v + gn0 − E

)(ε~k − ~~k · ~v + gn0 + E

)+ (gn0)

2 = 0(ε~k + gn0 +

(~~k · ~v − E

))(ε~k + gn0 −

(~~k · ~v − E

))− (gn0)2 = 0(

ε~k + gn0)2 − (~~k · ~v − E)2 − (gn0)