Plaintext
ΚΩΝΣΤΑΝΤΙΝΟΣ Ν. ΑΝΑΓΝΩΣΤΟΠΟΥΛΟΣ
Αναπληρωτής Καθηγητής
Εθνικό Μετσόβιο Πολυτεχνείο
Υπολογιστική Φυσική
Σημειώσεις
Το παρόν έργο αδειοδοτείται υπό τους όρους της άδειας Creative Commons Αναφορά
Δημιουργού - Μη Εμπορική Χρήση - Παρόμοια Διανομή 3.0. Για να δείτε ένα αντίγραφο της
άδειας αυτής επισκεφτείτε τον ιστότοπο
https://creativecommons.org/licenses/by-nc-sa/3.0/gr/
20 Απριλίου 2021
ii
ΠΕΡΙΕΧΟΜΕΝΑ
1 Κβαντική Δυναμική 1
1.1 Εισαγωγή . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Το Σχήμα του Visscher . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Το Πρόγραμμα . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4 Ελεύθερο Σωμάτιο . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.5 Σκέδαση . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.6 Ο Αρμονικός Ταλαντωτής . . . . . . . . . . . . . . . . . . . . . 22
1.7 Μετρήσεις . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.8 Παράρτημα . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
1.8.1 Απόδειξη Σχέσης (1.40) . . . . . . . . . . . . . . . . . . 32
1.8.2 Απόδειξη Σχέσης (1.48) . . . . . . . . . . . . . . . . . . 33
1.8.3 Αναμενόμενες Τιμές . . . . . . . . . . . . . . . . . . . . 36
1.9 Ασκήσεις . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
Bibliography 47
iii
iv ΠΕΡΙΕΧΟΜΕΝΑ
ΚΕΦΑΛΑΙΟ 1
Κβαντική Δυναμική
Στα προηγούμενα κεφάλαια συζητήθηκε ο αριθμητικός υπολογισμός του
ενεργειακού φάσματος και οι ιδιότητες των ιδιοκαταστάσεων της Χαμιλτονια-
νής ενός σωματιδίου που κινείται σε μία διάσταση. Οι καταστάσεις αυτές λέγο-
νται στάσιμες γιατί, όταν το σωματίδιο βρεθεί σε μία από αυτές, η κατάστασή
του δεν μεταβάλλεται με το χρόνο. Οποιαδήποτε άλλη κατάσταση μεταβάλλε-
ται δυναμικά με το χρόνο, σύμφωνα με την εξίσωση του Schrödinger. Η εξίσωση
αυτή είναι πρώτης τάξης ως προς το χρόνο και η γνώση της κατάστασης του
σωματιδίου μια χρονική στιγμή προσδιορίζει ντετερμινιστικά την κατάσταση
του σωματιδίου οποιαδήποτε άλλη χρονική στιγμή, υποθέτοντας ότι δεν γίνεται
μέτρηση. Οι στάσιμες καταστάσεις αποτελούν μία βάση στο χώρο Hilbert των
καταστάσεων του σωματιδίου και η χρονική εξέλιξη της κατάστασης του σωμα-
τιδίου μπορεί να εκφραστεί απλά ως γραμμικός συνδυασμός τους με τη χρονική
εξέλιξη να βρίσκεται μόνο στη χρονική εξάρτηση των συντελεστών του γραμ-
μικού συνδυασμού. Οι συντελεστές αυτοί μεταβάλλονται αρμονικά με το χρόνο
με συχνότητες που είναι ανάλογες της ενέργειας της αντίστοιχης ενεργειακής
ιδιοκατάστασης.
Στο κεφάλαιο αυτό θα παρουσιαστούν απλές αριθμητικές μέθοδοι για την
επίλυση της χρονοεξαρτώμενης εξίσωσης του Schrödinger. Η μέθοδος του Visscher
δίνει έναν αλγόριθμο ενός βήματος στο χρόνο, αλλά απαιτεί πολύ μικρό χρο-
νικό βήμα για να είναι ευσταθής. Η συνθήκη ευστάθειας είναι παρόμοια με τη
συνθήκη Courant (8.28) που συναντήσαμε στη μελέτη της εξίσωσης διάχυσης.
Στο επόμενο κεφάλαιο θα συζητηθεί το σχήμα Crank-Nicholson, το οποίο είναι
ευσταθές χωρίς να χρειάζεται να ικανοποιείται μια τέτοια συνθήκη, αλλά το
τίμημα είναι ότι σε κάθε βήμα χρειάζεται να αντιστρέφεται ένας τρι-διαγώνιος
πίνακας. Η απλότητα του σχήματος Visscher μας επιτρέπει να μελετήσουμε εύ-
κολα τη δυναμική συμπεριφορά ενός σωματιδίου περιορισμένου να κινείται στη
μία διάσταση και να καταλάβουμε βαθύτερα την κβαντική δυναμική.
1
2 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
1.1 Εισαγωγή
Η εξίσωση Schrödinger που δίνει τη χρονική εξέλιξη της κυματοσυνάρτησης
ψ(x, t) είναι
∂ψ(x, t)
ih̄ = Ĥψ(x, t) . (1.1)
∂t
Ĥ είναι ο τελεστής της Χαμιλτονιανής σε αναπαράσταση θέσης και δίνεται από
τη σχέση
h̄2 ∂ 2
Ĥ = − + V (x̂) . (1.2)
2m ∂x2
Η πυκνότητα πιθανότητας εύρεσης του σωματιδίου στη θέση x δίνεται από τη
σχέση
ρ(x, t) = |ψ(x, t)|2 = ψ(x, t)∗ ψ(x, t) , (1.3)
και η δυναμική εξέλιξη της ψ(x, t) σύμφωνα με την (1.1) εγγυάται ότι η συνολική
πιθανότητα ∫ +∞
P (t) = ψ(x, t)∗ ψ(x, t)dx (1.4)
−∞
παραμένει ανεξάρτητη του χρόνου, οπότε με κατάλληλη κανονικοποίηση της
ψ(x, t) μπορούμε να επιλέξουμε P (t) = P (0) = 1.
Από την (1.1), βλέπουμε ότι οι ιδιοκαταστάσεις της Χαμιλτονιανής
Ĥψn (x, t) = En ψn (x, t) , (1.5)
ικανοποιούν την εξίσωση
∂ψn (x, t)
ih̄ = En ψn (x, t) , (1.6)
∂t
οπότε η δυναμική τους εξέλιξη με το χρόνο δίνεται από τη σχέση
ψn (x, t) = ψn (x) e−iEn t/h̄ . (1.7)
Στις καταστάσεις αυτές, οι αναμενόμενες τιμές μιας οποιασδήποτε φυσικής πο-
σότητας που αναπαρίσταται από τον τελεστή Â είναι ανεξάρτητες του χρόνου:
∫ +∞
⟨Â⟩n = ψn (x, t)∗ A(x̂, p̂)ψn (x, t)dx
−∞
∫ +∞ (1.8)
∗ ∂
= ψn (x) A(x, −ih̄ )ψn (x)dx ,
−∞ ∂x
αφού (e−iEn t/h̄ )∗ e−iEn t/h̄ = 1.
1.1. ΕΙΣΑΓΩΓΗ 3
Οι ιδιοκαταστάσεις ψn (x) είναι μια βάση στον χώρο Hilbert των καταστά-
σεων που μπορεί να βρεθεί το σωματίδιο. Οπότε, αν τη χρονική στιγμή t = 0
το σωματίδιο βρίσκεται στην κατάσταση ψ(x, 0), αυτή μπορεί να γραφτεί σαν
γραμμικός συνδυασμός
∑
ψ(x, 0) = cn ψn (x) . (1.9)
n
Τότε η κατάσταση
∑
ψ(x, t) = cn ψn (x) e−iEn t/h̄ , (1.10)
n
ικανοποιεί την εξίσωση (1.1) με την αρχική συνθήκη (1.9), οπότε είναι η ζητού-
μενη λύση που δίνει τη δυναμική εξέλιξη της κατάστασης του σωματιδίου.
Η χρονική εξέλιξη της αναμενόμενης τιμής ⟨Â⟩(t) στην κατάσταση ψ(x, t)
δίνεται από τη σχέση
∫ +∞
⟨Â⟩(t) = ψ(x, t)∗ A(x̂, p̂)ψ(x, t)dx . (1.11)
−∞
Αντικαθιστώντας την (1.10) στην (1.11) παίρνουμε
∑
⟨Â⟩(t) = c∗n cm Ânm ei(En −Em )t/h̄ , (1.12)
n,m
όπου ∫ +∞
Ânm = ψn (x)∗ A(x̂, p̂) ψm (x)dx . (1.13)
−∞
Για να λύσουμε αριθμητικά την εξίσωση (1.1), είναι θεμιτό να την μετασχη-
ματίσουμε σε αδιάστατη μορφή. Για τον λόγο αυτό, επιλέγουμε μια κλίμακα μή-
κους L και εκφράζουμε όλες τις φυσικές ποσότητες σε μονάδες που προκύπτουν
από αυτή. Επιλέγουμε
h̄ p20 h̄2 h̄ mL2
p0 = , E0 = = , t0 = = , (1.14)
L m mL2 E0 h̄
ώς μονάδες ορμής, ενέργειας και χρόνου. Επίσης, είναι βολικό να θεωρήσουμε
τις μονάδες κυματαριθμού και ταχύτητας:
p0 1 p0 h̄
k0 =
= , v0 = = . (1.15)
h̄ L m mL
Μετρώντας τις φυσικές ποσότητες ως πολλαπλάσια των παραπάνω μονάδων,
μπορούμε να ορίσουμε τις αδιάστατες ποσότητες:
x p E t v
x̃ = , p̃ = , Ẽ = , t̃ = , ṽ = . (1.16)
L p0 E0 t0 v0
4 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
Η δυναμική ενέργεια V (x) μετριέται σε μονάδες ενέργειας E0 , οπότε ορί-
ζουμε την αδιάστατη συνάρτηση του δυναμικού
V (x̃L)
Ṽ (x̃) = (1.17)
E0
καθώς και την κυματοσυνάρτηση
√
ψ̃(x̃, t̃) = L ψ(x̃L, t̃ t0 ) . (1.18)
Οι κυματοσυναρτήσεις είναι κανονικοποιημένες στη μονάδα:
∫ +∞ ∫ +∞
∗
ψ(x, t) ψ(x, t)dx = ψ̃(x̃, t̃)∗ ψ̃(x̃, t̃)dx̃ = 1 . (1.19)
−∞ −∞
Αντικαθιστώντας τις (1.16), (1.17) και (1.18) στην (1.1), παίρνουμε:
∂ ψ̃(x̃, t̃) 1 ∂ 2 ψ̃(x̃, t̃)
i =− + Ṽ (x̃)ψ̃(x̃, t̃) . (1.20)
∂ t̃ 2 ∂ x̃2
Ας θεωρήσουμε το παράδειγμα του απλού αρμονικού ταλαντωτή:
1
V (x) = mω 2 x2 . (1.21)
2
√
Τότε, παίρνοντας ως L = h̄/mω, οι μονάδες μέτρησης είναι
√
h̄ √ 1
L= , p0 = h̄mω , E0 = h̄ω , t0 = , v 0 = p0 , (1.22)
mω ω
και το αδιάστατο δυναμικό
1
Ṽ (x̃) = x̃2 . (1.23)
2
Στις επόμενες παραγράφους θα χρησιμοποιούμε μόνο αδιάστατες ποσότη-
τες και θα παραλείπουμε τις περισπωμένες. Η εξίσωση (1.20) θα γράφεται για
συντομία:
∂ψ(x, t) 1 ∂ 2 ψ(x, t)
i =− + V (x)ψ(x, t) . (1.24)
∂t 2 ∂x2
Στα προβλήματα που θα λύσουμε, θα χρησιμοποιήσουμε ώς αρχική κατά-
σταση του σωματιδίου το κυματοπακέτο
1 (x−x0 )2
ik0 (x−x0 ) − 14
ψ(x, 0) = e e σ2 . (1.25)
(2πσ 2 )1/4
1.2. ΤΟ ΣΧΗΜΑ ΤΟΥ VISSCHER 5
Στην κατάσταση αυτή, οι δυνατές θέσεις παρατήρησης του σωματιδίου κατα-
νέμονται σύμφωνα με την πυκνότητα πιθανότητας
1 (x−x0 )2
− 21
ρ(x, 0) = ψ(x, 0)∗ ψ(x, 0) = e σ2 , (1.26)
(2πσ 2 )1/2
η οποία είναι Gaussian με κέντρο στο x0 και εύρος σ. Οι αναμενόμενες τιμές
⟨x⟩ = x0 , ⟨p⟩ = k0 , (1.27)
δείχνουν πως πρόκειται για σωματίδιο εντοπισμένο κοντά στη θέση x0 με μέση
ορμή k0 . Επίσης, είναι χρήσιμες και οι σχέσεις:
⟨x2 ⟩ = x20 + σ 2 , ⟨p2 ⟩ = k02 + 4σ1 2 ,
(1.28)
∆x2 = ⟨x2 ⟩ − ⟨x⟩2 = σ 2 , ∆p2 = ⟨p2 ⟩ − ⟨p⟩2 = 1
4σ 2
,
από όπου φαίνεται ότι η παράμετρος σ είναι ένα μέτρο της διασποράς των με-
τρήσεων της θέσης και της ορμής του σωματιδίου στην κατάσταση αυτή. Μικρό
σ σημαίνει μικρή διασπορά στη θέση και μεγάλη στην ορμή, και αντίστροφα.
1.2 Το Σχήμα του Visscher
Στην ενότητα αυτή παρουσιάζεται μια αριθμητική μέθοδος επίλυσης της
(1.1) χρησιμοποιώντας ένα σχήμα που δίνει την δυναμική εξέλιξη που προκύ-
πτει από την (1.1) με έναν αλγόριθμο ενός χρονικού βήματος (explicit method
ή leapfrog method). Ένας τέτοιος αλγόριθμος χρησιμοποιήθηκε στο Κεφάλαιο
8, ο οποίος συνοψίζεται από την εξίσωση (8.26). Η μέθοδος αυτή δεν μπορεί να
χρησιμοποιηθεί για τη λύση της (1.1) γιατί είναι ασταθής και η πιθανότητα (1.4)
δεν διατηρείται. Ο Visscher [1] πρότεινε μία βελτίωση της (8.26), όπου χρησιμο-
ποιεί το πραγματικό και το φανταστικό μέρος της κυματοσυνάρτησης ψ(x, t)
σε διαφορετικούς χρόνους. Αυτό έχει σαν αποτέλεσμα την ευστάθεια της δυνα-
μικής εξέλιξης και την διατήρηση της πιθανότητας (1.4).
Αναζητούμε αριθμητικές λύσεις της εξίσωσης (1.20)
∂ψ(x, t) 1 ∂ 2 ψ(x, t)
i =− + V (x)ψ(x, t) , (1.29)
∂t 2 ∂x2
στο διάστημα x ∈ [−xmax , xmax ] και t ∈ [0, tmax ]. Ορίζουμε το δισδιάστατο
πλέγμα, όπου σε μια συγκεκριμένη (διακριτή) χρονική στιγμή έχουμε Nx δια-
κριτά σημεία
2xmax
xi = −xmax + (i − 1)∆x , i = 1, . . . , Nx , ∆x = , (1.30)
Nx − 1
6 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
ενώ ο χρόνος παίρνει τις 2Nt διακριτές τιμές
tmax
tn = (n − 1)∆t , n = 1, . . . , Nt , ∆t = ,
Nt − 1
∆t
tn+ 1 = tn + , n = 1, . . . , Nt . (1.31)
2 2
Η κυματοσυνάρτηση ψ(x, t) χωρίζεται στο πραγματικό και το φανταστικό της
μέρος
ψ(x, t) = R(x, t) + i I(x, t) , (1.32)
όπου οι R(x, t) ≡ Re(ψ(x, t)) και I(x, t) ≡ Im(ψ(x, t)) είναι πραγματικές συ-
ναρτήσεις. Αντικαθιστώντας στην εξίσωση (1.85) παίρνουμε
∂R(x, t) 1 ∂ 2 I(x, t)
=− + V (x)I(x, t)
∂t 2 ∂x2
∂I(x, t) 1 ∂ 2 R(x, t)
=+ − V (x)R(x, t) (1.33)
∂t 2 ∂x2
①◆✍ ❂ ①✁✂✄
①✁✂✄ ✞ ✡①
①✁✂✄ ✞ ✷✡①
✡①
✌
✌
✌
✡t
✌
✌
✌
t✸☎✆ t✺☎✆ t✼☎✆
① ❂ ✵ ✿✿✿ ✿✿✿
t✶ ❂ ✵ t✆ t✸ t❥ ❂ ✭✝ ✞ ✟✠✡t t◆ ❂ t✁✂✄
✌
✌
✌
✌
✌
✌
✞①✁✂✄ ✰ ✷✡①
✞①✁✂✄ ✰ ✡①
①✶ ❂ ✞①♠☛☞
Σχήμα 1.1: Το χωροχρονικό πλέγμα στο σχήμα του Visscher αποτελείται από γεγονότα (xi , tn )
(μπλέ κύκλοι) και (xi , tn+ 21 ) (κόκκινοι κύκλοι), xi = −xmax + (i − 1)∆x, i = 1, . . . , Nx , tn =
(n−1)∆t, n = 1, . . . , Nt , tn+ 21 = tn +∆t/2. Στους μπλε κύκλους ορίζεται το πραγματικό μέρος
της κυματοσυνάρτησης Rin = R(xi , tn ), ενώ στους κόκκινους το φανταστικό Iin = I(xi , tn+ 21 ),
έτσι ώστε ψin = ψ(xi , tn ) ≡ Rin + iIin .
Η ιδέα του Visscher ήταν να θεωρήσει ένα σχήμα που να υπολογίζει τις
R(x, t) και I(x, t) σε διαφορετικές διαδοχικές χρονικές στιγμές tn και tn+ 1 . Το
2
1.2. ΤΟ ΣΧΗΜΑ ΤΟΥ VISSCHER 7
τμήμα του χωρόχρονου (x, t) με x ∈ [−xxmax , xxmax , ] και t ∈ [0, tmax ], προσεγγί-
ζεται από Nx × Nt διακριτά χωροχρονικά σημεία (xi , tn ) και (xi , tn+ 1 ), όπως
2
στο Σχήμα 1.1, έτσι ώστε
2xmax
xi = −xmax +(i − 1)∆x i = 1, . . . , Nx , ∆x = ,
Nx − 1
tmax
tn =(n − 1)∆t n = 1, . . . , Nt , ∆t = . (1.34)
Nt − 1
Στο όριο όπου ∆x → 0 και ∆t → 0, η R(x, t) θα προσεγγίζεται από τις
διακριτές τιμές {R(xi , tn )}, ενώ η I(x, t) από τις {I(xi , tn+ 1 )}. Τότε η κυματο-
2
συνάρτηση ψ(x, t) μπορεί να προσεγγιστεί από τις διακριτές τιμές¹
ψ(xi , tn ) ≡ R(xi , tn ) + i I(xi , tn+ 1 ) . (1.35)
2
Θα χρησιμοποιήσουμε τον συμβολισμό
ψin = ψ(xi , tn ) , Rin = R(xi , tn ) , Iin = I(xi , tn+ 1 ) , (1.36)
2
σύμφωνα με τον οποίο έχουμε
ψin = Rin + i Iin . (1.37)
Στο κεφάλαιο αυτό θα χρησιμοποιήσουμε τις συνοριακές συνθήκες
ψ(−xmax , t) = ψ(xmax , t) = 0 , (1.38)
οι οποίες, επί της ουσίας, βάζουν το σωματίδιο σε ένα απειρόβαθο πηγάδι δυ-
ναμικού
V (−xmax ) = V (xmax ) = +∞ . (1.39)
Τότε, το σχήμα Visscher δίνεται από τις επαγωγικές σχέσεις ενός χρονικού βή-
ματος
Rin+1 = Rin − a (Ii+1
n
− 2Iin + Ii−1
n
) + ∆t Vi Iin
Iin+1 = Iin + a (Ri+1
n+1
− 2Rin+1 + Ri−1 n+1
) − ∆t Vi Rin+1 , (1.40)
όπου
∆t
a= (1.41)
2∆x2
είναι η παράμετρος του Courant.
¹Για έναν πιο ακριβή ορισμό, δείτε την Άσκηση 3
8 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
Οι παραπάνω εξισώσεις λύνονται λαμβάνοντας υπόψη τις αρχικές συνθή-
κες και τις συνοριακές συνθήκες. Οι συνοριακές συνθήκες (1.38) υλοποιούνται
από τις σχέσεις
n
R1n = RN x
= I1n = INn x = 0 . (1.42)
Οι αρχικές συνθήκες καθορίζονται από την επιλεγμένη αρχική κυματοσυνάρ-
τηση ψ(x, 0). Λόγω της (1.35), η αρχική τιμή
ψi1 = Ri1 + i Ii1 , (1.43)
πρέπει να οριστεί με προσοχή. Σύμφωνα με τον ορισμό (1.36) Ii1 = I(xi , ∆t/2), η
οποία πρέπει να υπολογιστεί από την I(xi , 0) ≡ Im(ψ(xi , 0)). Κατ’ εξαίρεση, θα
χρησιμοποιήσουμε μία φορά την (1.40) για να προωθήσουμε την I(xi , 0) κατά
∆t/2, ορίζοντας:
a 1 ∆t
Ii1 = Im(ψ(xi , 0)) + (Ri+1 − 2Rn1 + Ri−1
1
)− Vi Ri1 , (1.44)
2 2
όπου Ri1 = Re(ψ(xi , 0)).
Η απόδειξη της εξίσωσης (1.40) μπορεί να βρεθεί στο Παράρτημα 1.8.1. Στην
εργασία [1] δείχνεται ότι η δυναμική εξέλιξη που δίνεται από τις εξισώσεις (1.40)
είναι ευσταθής αν
2 2 2
− ≤ Vi ≤ − (1.45)
∆t ∆t ∆x2
για όλα τα xi . Αυτό σημαίνει πως αν επιλέξουμε το ∆x, το βήμα ∆t θα πρέπει
να ικανοποιεί ταυτόχρονα τις παρακάτω ανισότητες:
2 2
∆t ≤ , ∆t ≤ , (1.46)
V+ + ∆x2 2 V−
όπου V+ = max{0, maxi {Vi }}, −V− = min{0, mini {Vi }}. Για V+ = 0, η αρι-
στερή ανισότητα δίνει τη συνθήκη a ≤ 1/2. Για V− = 0, η δεξιά ανισότητα δεν
περιορίζει το ∆t.
Όπως αναφέρθηκε και στην αρχή της παραγράφου, η πρόκληση για μια
αριθμητική μέθοδο που χρησιμοποιείται για την επίλυση της εξίσωσης Schrödinger
είναι να διατηρεί την ολική πιθανότητα (1.4) σταθερή στο χρόνο. Η πυκνότητα
πιθανότητας (1.3)
ρ(x, t) = ψ(x, t)∗ ψ(x, t) = R(x, t)2 + I(x, t)2 , (1.47)
δεν μπορεί να οριστεί ακριβώς στο σχήμα του Visscher, αφού οι συναρτήσεις
R(x, t) και I(x, t) ορίζονται στο πλέγμα σε διαφορετικές χρονικές στιγμές. Οι
ποσότητες
∑
Pn ≡ i ρni , ρni = (Rin )2 + Iin Iin−1
1 ∑ n+ 1
n+ 1 (1.48)
P n+ 2 ≡ i ρi 2 , ρi 2 = (Iin )2 + Rin Rin+1
1.3. ΤΟ ΠΡΟΓΡΑΜΜΑ 9
συγκλίνουν στις P (t) και ρ(t) στο όριο ∆t → 0 και διατηρούνται ακριβώς² από
τις εξισώσεις εξέλιξης (1.40) και τις συνοριακές συνθήκες (1.42). Το σφάλμα
στον υπολογισμό της P (t) στο όριο ∆t → 0, μετά από παρέλευση χρόνου tmax
είναι O(∆t).
1.3 Το Πρόγραμμα
Προκειμένου να προγραμματίσουμε τις επαγωγικές σχέσεις (1.40), αποθη-
κεύουμε στη μνήμη τις τιμές Rin , Rin+1 και Iin , Iin+1 του πραγματικού και φαντα-
στικού μέρους της κυματοσυνάρτησης τις χρονικές στιγμές tn , tn+1 και tn+ 1 ,
2
tn+1+ 1 αντίστοιχα. Για το λόγο αυτό ορίζουμε τα arrays
2
real ( 8 ) , d i m e n s i o n ( 2 , P ) : : repsi , impsi
στα οποία Rin → repsi(1,i), Rin+1 → repsi(2,i), Iin → impsi(1,i), Iin+1 →
impsi(2,i). Για ευκολία, ορίζουμε τις παραμέτρους
integer , parameter : : old = 1 , new =2
έτσι ώστε, λ.χ., Rin → repsi(old,i), Rin+1 → repsi(new,i).
Ο χρήστης δίνει στην είσοδο τις τιμές ∆x → dx, a = ∆t/2(∆x)2 → a,
xmax → xmax, tmax → tmax. Η τιμή του Nx = [2xmax /∆x] + 1, και θα πρέπει να
είναι μικρότερη από την επιλεγμένη τιμή της παραμέτρου P που καθορίζει τα
μεγέθη των arrays repsi και impsi. Αφού υπολογιστεί, η τιμή του ∆x επανα-
καθορίζεται, έτσι ώστε το 2xmax να είναι ακέραιο πολλαπλάσιό του. Η παράμε-
τρος ∆t → dt υπολογίζεται από τα ∆x, a, και στη συνέχεια δίνει τον αριθμό
των χρονικών βημάτων³ Nt = [tmax /∆t] + 1.
Η αρχική κυματοσυνάρτηση δίνεται από το κυματοπακέτο (1.25), και ο χρή-
στης δίνει στην είσοδο τις τιμές των παραμέτρων x0 → x0, σ → sigma και k0 →
k0, που είναι το κέντρο, η διασπορά και η μέση ορμή του κυματοπακέτου αντί-
στοιχα.
Οι τιμές του δυναμικού Vi → V(i) υπολογίζονται στο κυρίως πρόγραμμα
r e a l ( 8 ) , dimension ( P ) : : V
do ix = 1 , Nx
x = −xmax + ( ix −1) * dx
²Η απόδειξη δίνεται στο Παράρτημα 1.8.2.
³Επειδή το tmax μπορεί να μην είναι ακέραιο πολλαπλάσιο του ∆t, αφού υπολογιστεί το Nt ,
τα ∆t και a ξαναϋπολογίζονται από τις σχέσεις ∆t = tmax /(Nt − 1), a = ∆t/2(∆x)2 .
10 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
V ( ix ) = 0 . 5 D0 * x * x
end do
και κάθε δυναμικό πεδίο που μας ενδιαφέρει να μελετήσουμε θα πρέπει να προ-
γραμματιστεί ξεχωριστά.
Η δομή του προγράμματος είναι:
1. Είσοδος δεδομένων
2. Υπολογισμός δυναμικού
3. Υπολογισμός σταθερών παραμέτρων
4. Υπολογισμός αρχικής κυματοσυνάρτησης (1.25)
5. Υπολογισμός της Ii1 από την (1.44) (μισό βήμα ∆t/2)
6. Επαναληπτική εφαρμογή των σχέσεων (1.40) (από ένα βήμα ∆t)
7. Αποθήκευση των δεδομένων (n, tn , xi , ρni , Rin , Iin ) στο αρχείο psi.
Για εξοικονόμηση χώρου και ευκολότερη ανάλυση, η αποθήκευση γίνεται όταν
το n είναι πολλαπλάσιο μιας παραμέτρου tskip. Η παράμετρος tskip δίνεται
από τον χρήστη όταν ξεκινάει η εκτέλεση του προγράμματος.
Η τιμή της παραμέτρου του Courant a = ∆t/2(∆x)2 καθορίζεται ελεύθερα
από τον χρήστη, αλλά η subroutine compute_amax() υπολογίζει τη μέγιστη
τιμή amax που δίνει ευσταθείς λύσεις του σχήματος Visscher σύμφωνα με τις
(1.45). Το πρόγραμμα ελέγχει αν ισχύει a ≤ amax , και αν όχι, τυπώνει προειδο-
ποιητικό μήνημα στο stdout.
Το πλήρες πρόγραμμα για ένα σωμάτιο στο δυναμικό του αρμονικού ταλα-
ντωτή (1.23) δίνεται παρακάτω. Μπορεί να βρεθεί στο αρχείο tdse_fd.f90 του
συνοδευτικού λογισμικού.
! ============================================================
! Time Dependent S c h r o d i n g e r E q u a t i o n
!
! I m p l e m e n t s V i s s c h e r ’ s e x p l i c i t scheme
!
! f i l e : tdse_fd . f90
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
program tdse_fd
i m p l i c i t none
integer , parameter : : P = 100000
real ( 8 ) , d i m e n s i o n ( 2 , P ) : : repsi , impsi
real ( 8 ) , dimension ( P ) : : V , rho
1.3. ΤΟ ΠΡΟΓΡΑΜΜΑ 11
integer , parameter : : old = 1 , new =2
integer : : tskip
real (8) : : k0 , sigma , sigma2
real (8) : : xmax , tmax , x , t , dt , dx , a
real (8) : : amax , x0 , norm , pexp , delpsi
integer : : Nx , Nt
integer : : ix , it
real ( 8 ) , parameter : : tinypsi = 1 . 0 D−20
real ( 8 ) , parameter : : PI= a t a n 2 ( 1 . 0 D0 , 0 . 0 D0 ) * 2 . 0 D0
integer , parameter : : f_psi = 1 1 , f_rho =12
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
p r i n t * , ’ # E n t e r x0 , sigma , k0 , dx , a , xmax , tmax , t s k i p : ’
read * , x0 , sigma , k0 , dx , a , xmax , tmax , tskip
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
Nx = INT ( 2 . 0 D0 * xmax / dx ) + 1 ;
dx = 2 . 0 D0 * xmax / ( Nx − 1 )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! Potential :
V = 0 . 0 D0
do ix = 1 , Nx
x = −xmax + ( ix −1) * dx
V ( ix ) = 0 . 5 D0 * x * x
end do
c a l l compute_amax ( V , Nx , dx , amax ) ! Visscher ’ s conditions
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
dt = 2 . 0 D0 * a * dx * dx
Nt = INT ( tmax / dt ) + 1
dt = tmax / ( Nt −1)
a = 0 . 5 D0 * dt / ( dx * dx )
sigma2 = sigma * sigma
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
p r i n t * , ’ # −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− ’
print * , ’# x0 , sigma , k0 : ’
p r i n t * , ’ # ’ , x0 , sigma , k0
print * , ’# xmax , tmax , t s k i p : ’
p r i n t * , ’ # ’ , xmax , tmax , tskip
print * , ’# Nx , Nt , dx , d t : ’
p r i n t * , ’ # ’ , Nx , Nt , dx , dt
print * , ’# a , amax : ’
p r i n t * , ’ # ’ , a , amax
i f ( a>amax ) p r i n t * , ’ # WARNING : a >amax ’
p r i n t * , ’ # −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− ’
i f ( Nx > P ) s t o p ’ Nx > P ’
i f ( Nx < 3 ) s t o p ’ Nx t o o s m a l l ’
i f ( Nt / tskip < 2 ) s t o p ’ Nt t o o s m a l l ’
open ( u n i t =f_rho , f i l e = ’ rho ’ )
open ( u n i t =f_psi , f i l e = ’ p s i ’ )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! I n i t i a l condition :
12 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
t = 0 . 0 D0 ;
rho = 0 . 0 D0
repsi ( : , 1 ) = 0 . 0 D0 ; impsi ( : , 1 ) = 0 . 0 D0
repsi ( : , Nx ) = 0 . 0 D0 ; impsi ( : , Nx ) = 0 . 0 D0
norm = ( 2 . 0 D0 * PI * sigma2 ) * * ( − 0 . 2 5 D0 )
do ix = 2 , Nx−1
x = −xmax + ( ix −1) * dx
pexp = 0 . 2 5 D0 * ( x−x0 ) * ( x−x0 ) / sigma2
i f ( pexp < −l o g ( tinypsi ) ) t h e n ! a v o i d v e r y s m a l l numbers
repsi ( old , ix ) = norm * exp ( −pexp ) * c o s ( k0 * ( x−x0 ) )
impsi ( old , ix ) = norm * exp ( −pexp ) * s i n ( k0 * ( x−x0 ) )
else
repsi ( old , ix ) = 0 . 0 d0
impsi ( old , ix ) = 0 . 0 d0
end i f
rho ( ix ) = repsi ( old , ix ) * repsi ( old , ix ) &
+ impsi ( old , ix ) * impsi ( old , ix )
end do
norm = SUM ( rho ) * dx ! ~ ∫ rho ( x ) dx
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! Output i n i t i a l s t a t e :
it = 0 ; t= 0 . 0 D0
do ix = 1 , Nx
x = −xmax + ( ix −1) * dx
w r i t e ( f_psi , 1 0 1 ) it , t , x , rho ( ix ) ,&
repsi ( old , ix ) , impsi ( old , ix )
end do
w r i t e ( f_psi , * ) ’ ’ ;
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! Time e v o l u t i o n :
! F i r s t e v o l v e i m p s i by h a l f t i m e s t e p :
do ix = 2 , Nx−1
delpsi = a * ( repsi ( old , ix + 1 ) −2.0 D0 * repsi ( old , ix ) &
+ repsi ( old , ix −1) )−V ( ix ) * dt * repsi ( old , ix )
delpsi = 0 . 5 D0 * delpsi
i f ( a b s ( delpsi ) >tinypsi ) impsi ( old , ix ) =impsi ( old , ix ) + delpsi
end do
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! B e g i n l o o p on t i m e e v o l u t i o n :
do it = 1 , Nt−1
t = it * dt
! update r e p s i :
do ix = 2 , Nx−1
delpsi = −a * ( impsi ( old , ix + 1 ) −2.0 D0 * impsi ( old , ix ) &
+ impsi ( old , ix −1) ) +V ( ix ) * dt * impsi ( old , ix )
i f ( a b s ( delpsi ) >tinypsi ) repsi ( new , ix ) =repsi ( old , ix ) +delpsi
rho ( ix ) = repsi ( old , ix ) * repsi ( new , ix ) &
+ impsi ( old , ix ) * impsi ( old , ix )
end do
1.3. ΤΟ ΠΡΟΓΡΑΜΜΑ 13
do ix = 2 , Nx−1
delpsi = a * ( repsi ( new , ix + 1 ) −2.0 D0 * repsi ( new , ix ) &
+ repsi ( new , ix −1) )−V ( ix ) * dt * repsi ( new , ix )
i f ( a b s ( delpsi ) >tinypsi ) impsi ( new , ix ) =impsi ( old , ix ) +delpsi
end do
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! p r i n t norm2 o f p s i : ∫ rho ( x ) dx
norm = SUM ( rho ) * dx
w r i t e ( f_rho , * ) it , t , norm
! print psi :
i f (MOD( it , tskip ) == 0 ) t h e n
do ix = 1 , Nx
x = −xmax + ( ix −1) * dx
w r i t e ( f_psi , 1 0 1 ) it , t , x , rho ( ix ) , repsi ( new , ix ) , impsi ( new , ix )
end do
w r i t e ( f_psi , * ) ’ ’
end i f
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! i n t e r c h a n g e new<−> o l d v a l u e s :
repsi ( old , : ) = repsi ( new , : )
impsi ( old , : ) = impsi ( new , : )
end do ! do i t = 1 , Nt end t i m e e v o l u t i o n
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
c l o s e ( f_rho )
c l o s e ( f_psi )
1 0 1 FORMAT ( I20 , 1 0 0 0 E30 . 1 5 )
end program tdse_fd
! ============================================================
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
subroutine compute_amax ( V , Nx , dx , a )
i m p l i c i t none
integer : : Nx
r e a l ( 8 ) , d i m e n s i o n ( Nx ) : : V
real (8) : : dx , a
real (8) : : V1 , V2 , dt1 , dt2 , dt
V1 = MAXVAL ( V ( 1 : Nx ) ) ; i f ( V1 <= 0 . 0 D0 ) V1 = 0 . 0 D0
V2 = −MINVAL ( V ( 1 : Nx ) ) ; i f ( V2 <= 0 . 0 D0 ) V2 = 0 . 0 D0
dt1 = 2 . 0 D0 / ( V1 + ( 2 . 0 D0 / ( dx * dx ) ) )
dt2 = HUGE ( dt )
i f ( V2 > 0 . 0 D0 ) dt2 = 2 . 0 D0 / V2
dt = dt2
i f ( dt1 <dt2 ) dt = dt1
a = 0 . 5 D0 * dt / ( dx * dx )
14 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
end s u b r o u t i n e compute_amax
1.4 Ελεύθερο Σωμάτιο
Στην παράγραφο αυτή θα μελετήσουμε το σωμάτιο στο δυναμικό
V (x) = 0 . (1.49)
Στην πραγματικότητα, λόγω των συνοριακών συνθηκών (1.38), το σωμάτιο βρί-
σκεται σε ένα απειρόβαθο πηγάδι δυναμικού (1.39) πλάτους 2xmax , αλλά όσο το
σωμάτιο είναι βρίσκεται μακριά από τα τοιχώματα (δηλ. η τιμή της ψ(x, t) είναι
αμελητέα στη γειτονιά των σημείων x = ±xmax ), η δυναμική του σωματιδίου θα
προσεγγίζει πολύ καλά την δυναμική του ελεύθερου σωματιδίου.
Αφού θέσουμε V(ix) = 0.0D0 στο πρόγραμμα που βρίσκεται στο αρχείο
tdse_fd.f90, μεταγλωττίζουμε και εκτελούμε το πρόγραμμα με τις εντολές:
> g f o r t r a n tdse_fd . f90 −o t
> echo −7 0.5 3 . 0 0 . 0 2 0 . 4 5 10 2.5 200 |./t
# E n t e r x0 , sigma , k0 , dx , a , xmax , tmax , t s k i p :
# −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
# x0 , sigma , k0 :
# −7.0 0 . 5 3.0
# xmax , tmax , t s k i p :
# 10.0 2 . 5 0 200
# Nx , Nt , dx , dt :
# 1 0 0 1 6 9 4 5 2 . 0 E−002 3 . 6 0 0 2 3 E−004
# a, amax :
# 0.450029 0.50
# −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
Στο αρχείο rho βρίσκουμε σε τρεις στήλες τις τιμές (n, tn , P n ). Παρατηρήστε
ότι οι τιμές P n παραμένουν σταθερές⁴.
Στο αρχείο psi βρίσκουμε σε 6 στήλες τις τιμές (n, tn , xi , ρni , Rin , Iin ), Κάθε
φορά που αλλάζει η τιμή του n έχουμε τυπώσει μια κενή γραμμή. Έτσι, μπο-
ρούμε με την awk να βρούμε πόσα στιγμιότυπα της κυματοσυνάρτησης έχουμε
συλλέξει:
> awk ’ NF < 1 { n + + }END { p r i n t ” Number o f f r a m e s = ” , n } ’ psi
Number of frames= 35
⁴Για άλλες επιλογές παραμέτρων, λ.χ. μεγαλύτερο tmax, μπορεί, λόγω αριθμητικών σφαλ-
μάτων, να υπάρχει μια μικρή μεταβολή στα τελευταία σημαντικά ψηφία του P n
1.4. ΕΛΕΥΘΕΡΟ ΣΩΜΑΤΙΟ 15
Θυμηθείτε πως η εσωτερική μεταβλητή NF της awk μετράει τον αριθμό των λέ-
ξεων (NF = number of fields) που έχει κάθε γραμμή του αρχείου που διαβάζει.
Στην περίπτωση του αρχείου psi, οι γραμμές που περιέχουν δεδομένα έχουν έξι
“λέξεις” (NF=6), ενώ σε κάθε αλλαγή στιγμιότυπου αφήνουμε μια κενή γραμμή
με μηδέν “λέξεις” (NF=0).
Φυσικά, ο λόγος που ο αριθμός των στιγμιότυπων είναι μικρότερος του Nt =
6945, είναι ότι τα τυπώνουμε κάθε φορά που το n είναι πολλαπλάσιο της πα-
ραμέτρου tskip, την οποία στο παραπάνω παράδειγμα έχουμε θέσει να είναι
200, οπότε τυπώνονται για nframe = 0, 1, . . . , [Nt /ttmax ] = 0, 1, . . . , 34.
Με τον παραπάνω τρόπο μπορούμε να εκτυπώσουμε το στιγμιότυπο της κυ-
ματοσυνάρτησης που βρίσκεται στο αρχείο psi για tn αμέσως μετά από μια
επιλεγμένη τιμή του χρόνου t:
> awk −v t =2 ’ $2 >t { p = 1 ; p r i n t } p==1&&NF < 1 { e x i t } ’ psi > psi_t2 . 0
Με την εντολή αυτή, στο αρχείο psi_t2.0 θα βρούμε το στιγμιότυπο ψin για
tn ≈ 2.016, το οποίο είναι το αμέσως επόμενο στιγμιότυπο του t = 2.
Για να δούμε γραφικά τη χρονική εξέλιξη της ρ(x, t) μπορούμε να κάνουμε
την τρισδιάστατη γραφική παράσταση:
gnuplot > splot ” p s i ” using 2 : 3 : 4 with dots
t=0.0
t=0.3
0.7
t=1.0
0.8 t=2.4
0.7 0.8 0.6
0.6 0.7
0.5 0.6
0.4 0.5 0.5
0.3 0.4
0.2 0.3
ρ(x,t)
0.1 0.2 0.4
0 0.1
-0.1 0
-0.1 0.3
0 0.2
0.5
1
1.5
2 0 5 10 0.1
2.5 -10 -5
0
-10 -5 0 5 10
x
Σχήμα 1.2: Η πυκνότητα πιθανότητας ρ(x, t) για το ελεύθερο σωμάτιο, όπως υπολογίστηκε από
το πρόγραμμα tdse_fd.f90 στη σελίδα 14.
Το αποτέλεσμα θα είναι⁵ το αριστερό σχήμα του Σχήματος 1.2
Πιο απλό είναι να κάνουμε τη γραφική παράσταση ενός στιγμιότυπου της
κυματοσυνάρτησης. Αυτό μπορεί να γίνει αποθηκεύοντας στιγμιότυπα σε ξεχω-
ριστά αρχεία και κάνοντας τη γραφική τους παράσταση με την awk, όπως δεί-
⁵Για να γίνει όπως στο Σχήμα 1.2 έχουμε κάνει splot with lines αφού δώσαμε τις εντολές
set hidden3d;set pm3d;, αλλά αυτό θα σας απορροφήσει πολλούς υπολογιστικούς πόρους.
16 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
ξαμε παραπάνω, ή επιλέγοντας τα στιγμιότυπα με το κατάλληλο φίλτρο μέσα
από το gnuplot:
gnuplot > p l o t \
” <awk ’ $2 > 0 . 0 { p = 1 ; print } p==1&&NF < 1 { exit } ’ psi ” u 3:4 w l,\
” <awk ’ $2 > 0 . 3 { p = 1 ; print } p==1&&NF < 1 { exit } ’ psi ” u 3:4 w l,\
” <awk ’ $2 > 1 . 0 { p = 1 ; print } p==1&&NF < 1 { exit } ’ psi ” u 3:4 w l,\
” <awk ’ $2 > 2 . 4 { p = 1 ; print } p==1&&NF < 1 { exit } ’ psi ” u 3:4 w l
Η παραπάνω εντολή θα φτιάξει τη δεξιά γραφική παράσταση του Σχήματος
1.2.
0.4 2.5
ρ(x,t)
Re(ψ(x,t))
0.3
Im(ψ(x,t))
0.2 2
0.1
Δx(t)
1.5
0
-0.1
1
-0.2
-0.3
0.5
-0.4
-10 -5 0 5 10 0 0.5 1 1.5 2 2.5
x t
Σχήμα 1.3: Αριστερά φαίνονται οι ρ(x, t), R(x, t), I(x, t), για t = 2.4 για το ελεύθερο σωμάτιο,
όπως υπολογίστηκε από το πρόγραμμα tdse_fd.f90 στη σελίδα 14. Δεξιά υπολογίζεται η δια-
σπορά ∆x(t) (διακριτά σημεία) τα οποία συγκρίνονται με τη σχέση (1.50) και την ασυμπτωτική
ευθεία t/(2σ).
Με παρόμοιες εντολές μπορούμε να φτιάξουμε και τις γραφικές παραστά-
σεις του πραγματικού και του φανταστικού μέρους της κυματοσυνάρτησης για
κάποια χρονική στιγμή. Με τις παρακάτω εντολές μπορούμε να φτιάξουμε τη
γραφική παράσταση του Σχήματος (1.3)
gnuplot > p l o t \
” <awk ’ $2 > 2 . 4 { p = 1 ; p r i n t } p==1&&NF < 1 { e x i t } ’ p s i ” u 3 : 4 w l , \
” <awk ’ $2 > 2 . 4 { p = 1 ; p r i n t } p==1&&NF < 1 { e x i t } ’ p s i ” u 3 : 5 w l , \
” <awk ’ $2 > 2 . 4 { p = 1 ; p r i n t } p==1&&NF < 1 { e x i t } ’ p s i ” u 3 : 6 w l
Με τους παραπάνω υπολογισμούς μπορούμε να δούμε μια βασική ιδιότητα
της κατάστασης ενός ελεύθερου σωματιδίου που δίνεται από ένα εντοπισμένο
κυματοπακέτο, όπως αυτό της εξίσωσης (1.25). Το κυματοπακέτο εξαπλώνε-
ται απεριόριστα με το χρόνο και το σωματίδιο εντοπίζεται σε ολοένα και με-
γαλύτερο διάστημα. Αυτό εκφράζεται ποσοτικά από τη διασπορά ∆x(t) =
1.4. ΕΛΕΥΘΕΡΟ ΣΩΜΑΤΙΟ 17
√
⟨x2 ⟩(t) − (⟨x⟩(t))2 , η οποία μπορεί να δειχθεί ότι είναι⁶
√
t 4σ 4
∆x(t) = 1+ 2 . (1.50)
2σ t
Για αρκετά μεγάλους χρόνους, η διασπορά αυξάνεται ανάλογα με τον χρόνο
και αυτό μπορούμε να το δούμε στην αριστερή γραφική παράσταση του Σχή-
ματος 1.4 όπου φαίνεται η εξάπλωση του κυματοπακέτου. Η σχέση μπορεί να
ελεχθεί και ποσοτικά αν κάνουμε έναν γρήγορο υπολογισμό⁷∫ της αβεβαιότη-
+xmax
τας ∆x(t). Από τη σχέση (1.11) έχουμε τις σχέσεις ⟨x⟩ = −xmax xρ(x, t)dx,
∫ +x
⟨x2 ⟩ = −xmax x2 ρ(x, t)dx, και το παρακάτω πρόγραμμα, γραμμένο σε awk⁸, υπο-
max
λογίζει τις τιμές (t, ⟨x⟩(t), ⟨x2 ⟩(t), ∆x(t)):
# ! / u s r / b i n / awk −f
# −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
# When t h e l i n e h a s d a t a : NF>1 ( i n our c a s e NF = 6 )
NF > 1 {
t =$2 ; # current time
x +=$3 * $4 ; # Σ_i x * ρ (x , t )
x2+=$3 * $3 * $4 ; # Σ _ i x * x * ρ ( x , t )
}
# A blank l i n e s e p a r a t e s d a t a s e t s of d i f f e r e n t time :
NF < 1 {
dx = 0 . 0 2 ; # p u t h e r e t h e v a l u e o f dx
x *= dx ; # ~ ∫ x * ρ ( x , t ) * dx
x2 *= dx ; # ~ ∫ x * x * ρ ( x , t ) * dx
DX = s q r t ( x2−x * x ) ; # Δx = s q r t ( <x * x> − <x >* < x> )
p r i n t t , x , x2 , DX ;
x = 0 ; x2 = 0 ; # r e s e t sums f o r t h e n e x t d a t a s e t
}
Το παραπάνω πρόγραμμα μπορείτε να το βρείτε στο συνοδευτικό λογισμικό στο
αρχείο tdse_fd_DX.awk. Αν το γράψετε σε δικό σας αρχείο, θα πρέπει πρώτα
να του δώσετε άδεια εκτέλεσης:
chmod a+x tdse_fd_DX . awk
και μετά να το τρέξετε:
⁶Δείτε τις εξισώσεις (1.93).
⁷Σε επόμενη ενότητα θα κάνουμε πιο συστηματική μελέτη υπολογισμού αναμενόμενων τι-
μών.
⁸Το πρόγραμμα αυτό μπορεί να γραφτεί σε μια γραμμή awk 'NF>1 t=$2;
x+=$3*$4; x2+=$3*$3*$4 NF<1 dx=0.02; x*=dx; x2*=dx; DX=sqrt(x2-x*x) ;print
t,x,x2,DX; x=0; x2=0' και να χρησιμοποιηθεί από τη γραμμή εντολών. Αν προτιμάτε,
μπορείτε να ενσωματώσετε τον υπολογισμό στον πηγαίο κώδικα tdse_fd.f90.
18 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
> . / tdse_fd_DX . awk psi > DX . dat
gnuplot > plot ”DX . d a t ” using 1 : 4 with lines
Το αποτέλεσμα φαίνεται στην δεξιά γραφική παράσταση του Σχήματος 1.3. Πα-
ρόλη την απλοϊκότητα του υπολογισμού, η συμφωνία με τη σχέση (1.50) είναι
πολύ καλή. Στο ίδιο σχήμα φαίνεται η ασυμπτωτική συμπεριφορά ∆x(t) ∼
t/(2σ) καθώς το t → ∞. Για να φτιάξετε το σχήμα αυτό, χρησιμοποιώντας
το πρόγραμμα που γράψατε στο αρχείο tdse_fd_DX.awk, εκτελέστε τις παρα-
κάτω εντολές στο gnuplot:
gnuplot > DT ( t ) = s q r t ( s * s + ( t / ( 2 * s ) ) * * 2 )
gnuplot > asympt ( t ) = t / ( 2 * s )
gnuplot > s = 0.5
gnuplot > p l o t ” < . / t d s e _ f d _ D X . awk p s i ” u 1 : 4 , DT ( x ) , asympt ( x )
0.07
t=0.00 t=10.9
t=0.43 t=14.0
0.7 0.06
t=0.86 t=29.9
t=1.94
0.6 t=4.03
0.05
t=6.60
0.5 t=8.90
0.04
ρ(x,t)
ρ(x,t)
0.4
0.03
0.3
0.02
0.2
0.1 0.01
0 0
-30 -20 -10 0 10 20 30 -30 -20 -10 0 10 20 30
x x
Σχήμα 1.4: Η πυκνότητα πιθανότητας ρ(x, t), όπως υπολογίστηκε από το πρόγραμμα
tdse_fd.f90 για Vi = 0. Στο αριστερό σχήμα βλέπουμε την εξάπλωση του κυματοπακέτου
όσο το σωμάτιο κινείται, πρακτικά ελεύθερο, καθώς περνάει ο χρόνος. Στο δεξί σχήμα βλέπουμε
το αποτέλεσμα των ανακλάσεων της κυματοσυνάρτησης στα τοιχώματα V (±xmax ) = +∞,
xmax = 30.
Όπως είπαμε και παραπάνω, λόγω των συνοριακών συνθηκών (1.38), το σω-
μάτιο δεν είναι ελεύθερο, αλλά βρίσκεται σε ένα απειρόβαθο πηγάδι δυναμικού
(1.39) πλάτους 2xmax . Άρα από κάποια χρονική στιγμή και μετά⁹ η κυματοσυ-
νάρτηση θα ανακλάται στα τοιχώματα και τα κύματα που κινούνται σε αντίθε-
τες κατευθύνσεις θα συμβάλλουν. Το αποτέλεσμα της συμβολής φαίνεται στη
δεξιά γραφική παράσταση του Σχήματος 1.4.
⁹Εκτιμήστε την από το γεγονός ότι, όσο το σωμάτιο είναι (κατά προσέγγιση) ελεύθερο,
⟨v⟩(t) = ⟨p⟩(t) = k0 .
1.5. ΣΚΕΔΑΣΗ 19
1.5 Σκέδαση
Στην ενότητα αυτή θα μελετήσουμε τη σκέδαση ενός σωματιδίου σε ένα τε-
τραγωνικό δυναμικό. Θα ξεκινήσουμε με τη μελέτη ενός φράγματος που δίνεται
από το δυναμικό: {
V0 |x| < b
V (x) = , (1.51)
0 |x| > b
με V0 > 0. Όπως γνωρίζουμε από την κβαντομηχανική, μέρος της κυματοσυ-
νάρτησης ανακλάται και μέρος της συνεχίζει να διαδίδεται. Στην περίπτωση
που το σωμάτιο έχει καθορισμένη ενέργεια E, μη μηδενική διάδοση συμβαίνει
ακόμα και αν E < V0 λόγω του φαινομένου σήραγγας, ενώ μη μηδενική ανά-
κλαση συμβαίνει ακόμα και αν E > V0 . Το κυματοπακέτο που δίνεται από τη
σχέση (1.25) είναι υπέρθεση ιδιοκαταστάσεων όλων των δυνατών ενεργειών¹⁰,
οπότε αναμένεται μέρος του να συνεχίσει τη διάδοση και μέρος του να ανακλα-
στεί.
Θα μελετήσουμε την περίπτωση V0 = 6, b = 0.5. Μεταβάλλουμε τον κώδικα
στο αρχείο tdse_fd.f90 στο σημείο που υπολογίζεται το δυναμικό ως εξής:
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! Potential :
V = 0 . 0 D0
do ix = 1 , Nx
x = −xmax + ( ix −1) * dx
i f ( ABS ( x ) <= 0 . 5 D0 ) V ( ix ) = 6 . 0 D0
end do
Για να δούμε κυρίως τα φαινόμενα σκέδασης από το φράγμα (1.51), πρέ-
πει να επιλέξουμε τις παραμέτρους x0 , σ, xmax , tmax , έτσι ώστε να ελαχιστοποι-
ήσουμε τα φαινόμενα της αλληλεπίδρασης (ανάκλαση/συμβολή) του κυματο-
πακέτου με τα τοιχώματα του πηγαδιού στα x = ±xmax . Επιλέγοντας
> echo −6 0 . 5 3.0 0.02 0.45 300 40 1000 |./t
# E n t e r x0 , sigma , k0 , dx , a , xmax , tmax , t s k i p :
....
μπορούμε να δούμε τη σκέδαση του αρχικού κυματοπακέτου, όπως φαίνεται
στα στιγμιότυπα του Σχήματος 1.5. Με την επιλογή k0 = 3, η μέση αρχική ορμή
του σωματιδίου ⟨p⟩ = ⟨v⟩ = 3, η μέση ορμή του τετραγώνου της ορμής ⟨p2 ⟩ =
k02 + 1/(4σ 2 ) = 10, άρα η μέση τιμή της ενέργειας είναι ⟨H⟩ = ⟨p2 ⟩/2 = 5.
¹⁰Φυσικά δεν συνεισφέρουν όλες οι ενεργειακές ιδιοκαταστάσεις το ίδιο, κυρίως συνεισφέ-
ρουν αυτές με ενεργεια γύρω από τη μέση τιμή ⟨H⟩ = ⟨p2 ⟩/2 που δίνεται από την (1.28).
20 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
0.8 t=0.0 t=2.5
0.6
V(x)/10 V(x)/10
0.7
0.5
0.6
0.4
0.5
ρ(x,t)
ρ(x,t)
0.4 0.3
0.3
0.2
0.2
0.1
0.1
0 0
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
x x
t=3.5 t=18.5
0.6 0.02
V(x)/10 t=39.0
0.5
0.015
0.4
ρ(x,t)
ρ(x,t)
0.3 0.01
0.2
0.005
0.1
0 0
-15 -10 -5 0 5 10 15 -200 -100 0 100 200
x x
Σχήμα 1.5: Σκέδαση από το φράγμα δυναμικού που δίνεται από την εξίσωση (1.51) με V0 =
6, b = 0.5. Επάνω αριστερά φαίνεται η αρχική κατάσταση πριν την σκέδαση. Τα δύο επόμενα
σχήματα δείχνουν δύο στιγμιότυπα της σκέδασης, ενώ κάτω δεξιά τα διαδιδόμενα και ανακλό-
μενα κυματοπακέτα έχουν πρακτικά χωριστεί.
Άρα ⟨H⟩ < V0 = 6, τιμή η οποία θα μειωθεί περαιτέρω με το άπλωμα του
κυματοπακέτου (αύξηση σ). Αν θεωρήσουμε ένα ανάπτυγμα του κυματοπακέτου
στη βάση των ενεργειακών ιδιοκαταστάσεων, τότε οι συνιστώσες με ενέργεια
E < V0 θα συνεισφέρουν στη διάδοση μέσω του φαινομένου της σήραγγας
(εκθετική μείωση). Οι υπόλοιπες,
√ με E > V0 , θα διαδίδονται μέσα στο φράγμα
με μήκος κύματος k = 2π/λ = 2(E − V0 ).
Μετά από αρκετό χρόνο από τη σκέδαση, τα ανακλώμενα και διαδιδόμενα
κυματοπακέτα διαδίδονται πρακτικά ανεξάρτητα το ένα από το άλλο. Τότε
μπορούμε να ρωτήσουμε ποια είναι η πιθανότητα το σωματίδιο να έχει ανακλα-
στεί ή να έχει περάσει το φράγμα. Αν pR,∞ είναι η πιθανότητα να ανακλαστεί
και pT,∞ να το περάσει, τότε θα έχουμε
∫ 0
pR,∞ = lim ρ(x, t)dx
t→+∞ −x
∫ +xmaxmax
pT,∞ = lim ρ(x, t)dx . (1.52)
t→+∞ 0
Το παρακάτω πρόγραμμα, γραμμένο
∑ ∑ σε awk, υπολογίζει τα αθροίσματα pR (tn ) ≈
xi <0 ρi ∆x και pT (tn ) ≈ xi >0 ρi ∆x, τα οποία μπορούμε να δούμε γραφικά,
n n
όπως στο Σχήμα 1.6:
1.5. ΣΚΕΔΑΣΗ 21
1 1
pR pR,∞
pT pT,∞
pR,∞
0.8 pT,∞ 0.8
0.6 0.6
0.4
0.4
0.2
0.2
0
0 0 2 4 6 8 10 12 14 16
0 1 2 3 4 5 6
V0
t
Σχήμα 1.6: (Αριστερά) Σκέδαση από το φράγμα δυναμικού που δίνεται από την εξίσωση (1.51)
∫0 ∫ +x
με V0 = 6, b = 0.5, k0 = 3, σ = 0.5. pR (t) = −xmax ρ(x, t)dx, pT (t) = 0 max ρ(x, t)dx,
pR,∞ = pR (+∞) = 0.743346, pT,∞ = pT (+∞) = 0.256654. (Δεξιά) Η εξάρτηση των pR,∞ ,
pT,∞ από το ύψος του φράγματος V0 . Όλες οι άλλες παράμετροι είναι οι ίδιες.
# ! / u s r / b i n / awk −f
# −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
NF > 1 {
t =$2 ;
x =$3 ;
rho=$4 ;
i f ( x < 0 ) R += rho ;
i f ( x > 0 ) T += rho ;
}
# A blank l i n e s e p a r a t e s d a t a s e t s of d i f f e r e n t time :
NF < 1 {
dx = 0 . 0 2 ; # p u t h e r e t h e v a l u e o f dx
R *= dx ; # ~ ∫ρ ( x , t ) * dx
T *= dx ;
print t , R , T ;
R=0;T=0; # r e s e t sums f o r t h e n e x t d a t a s e t
}
Στο Σχήμα αυτό βλέπουμε πως τα pR,T (tn ) συγκλίνουν για t > 2.2 και μπο-
ρούμε εύκολα να υπολογίσουμε τα όρια pR,∞ και pT,∞ . Συμπεραίνουμε πως η
πιθανότητα το σωμάτιο να περάσει το φράγμα μετά από αρκετά μεγάλο χρόνο
είναι περίπου 1/4.
Για δεδομένο αρχικό κυματοπακέτο, οι πιθανότητες pR,∞ και pT,∞ εξαρτώ-
νται από τις παραμέτρους V0 και b. Κρατώντας το b σταθερό, η αύξηση του V0
οδηγεί σε αύξηση του pR,∞ , μέχρι να έχουμε σχεδόν ολική ανάκλαση. Αυτό μπο-
ρεί να γίνει κατανοητό αν θεωρήσουμε το ανάπτυγμα της κυματοσυνάρτησης
σε γραμμικό συνδυασμό ενεργειακών ιδιοκαταστάσεων, όπως στην εξίσωση
(1.10). Οι καταστάσεις με ενέργεια E < V0 διαδίδονται μέσω του φράγματος
λόγω του φαινομένου σήραγγας. Οι κυματοσυναρτήσεις τους έχουν μέτρο που
22 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
−kd
ελαττώνεται
√ εκθετικά καθώς διειδύει σε απόσταση d μέσα στο φράγμα ως e
με k = 2(V0 − E). Άρα, για τις καταστάσεις αυτές, το μήκος διείσδυσης είναι
dmax ∼ 1/k, το οποίο μειώνεται καθώς αυξάνει το V0 . Επίσης, όταν ⟨H⟩ ≪ V0 ,
οι καταστάσεις με E > V0 έχουν αμελητέες συνιστώσες στην εξίσωση (1.10)
και το φαινόμενο σήραγγας είναι αυτό που επικρατεί στη διάδοση. Αντίθετα,
μείωση του V0 οδηγεί σε αύξηση του pT,∞ , έως να έχουμε V0 = 0, οπότε έχουμε
ολική διάδοση. Στο δεξί Σχήμα 1.6 βλέπουμε την εξάρτηση αυτή από το V0 .
Με τον ίδιο τρόπο, κρατώντας το V0 σταθερό, η αύξηση του b μειώνει την
pT,∞ , μέχρι οι συνιστώσες με E < V0 σχεδόν να μην διαδίδονται, και να έχουμε
διάδοση μόνο των E > V0 συνιστωσών.
0.035
0.8 t=0.00 t=10
t=0.72 0.03 t=28
t=1.44 V(x)/1000
t=2.52 0.025
0.6 t=3.60
V(x)/100 0.02
0.015
0.4
ρ(x,t)
ρ(x,t)
0.01
0.2 0.005
0
0 -0.005
-0.01
-0.2 -0.015
-15 -10 -5 0 5 10 15 -150 -100 -50 0 50 100 150
x x
Σχήμα 1.7: Σκέδαση από το πηγάδι δυναμικού που δίνεται από την εξίσωση (1.51) με V0 =
−12, b = 2.
Τέλος, μπορούμε να μελετήσουμε και τη σκέδαση από ένα πηγάδι δυναμι-
κού, θέτοντας στην (1.51) V0 < 0. Στο Σχήμα 1.7 φαίνονται στιγμιότυπα της
πυκνότητα πιθανότητας ρ(x, t), όταν V0 = −12, b = 2. Αν θεωρήσουμε ένα
ανάπτυγμα (1.10) του κυματοπακέτου στη βάση των ενεργειακών ιδιοκαταστά-
σεων, τότε μέσα στο πηγάδι√οι συνιστώσες με ενέργεια E θα διαδίδονται με
μήκος κύματος k = 2π/λ = 2(E − V0 ).
1.6 Ο Αρμονικός Ταλαντωτής
Ο απλός αρμονικός ταλαντωτής, όταν τα μεγέθη κανονικοποιηθούν στις
μονάδες που δίνονται από τις (1.22), δίνεται από τη συνάρτηση δυναμικού
1
V (x) = x2 . (1.53)
2
Αφού προγραμματίσουμε τις τιμές V(ix) στο αρχείο tdse_fd.f90, όπως δεί-
ξαμε στη σελίδα 9, και μεταγλωττίσουμε, μπορούμε εκτελέσουμε την εντολή:
> echo 0 0.5 3 0 . 0 2 0 . 4 5 15 7 40 |./t
1.6. Ο ΑΡΜΟΝΙΚΟΣ ΤΑΛΑΝΤΩΤΗΣ 23
# E n t e r x0 , sigma , k0 , dx , a, xmax , tmax , t s k i p :
...
Στη συνέχεια μπορούμε να μελετήσουμε γραφικά την κυματοσυνάρτηση και την
πυκνότητα πιθανότητας, όπως και στις προηγούμενες ενότητες. Παρατηρούμε
ότι η ρ(x, t) είναι περιοδική συνάρτηση του χρόνου, με περίοδο που καθορίζεται
από την επιλογή ω = 1. Η διασπορά του κυματοπακέτου δεν αυξάνει διαρκώς με
τον χρόνο, όπως στην περίπτωση του ελεύθερου σωμάτιου, αλλά είναι και αυτή
περιοδική συνάρτηση του χρόνου. Είναι ελάχιστη στο σημείο ισορροπίας και
μέγιστη στη θέση μέγιστης απομάκρυνσης. Από την αρχή απροσδιοριστίας του
Heisenberg συμπεραίνουμε πως η διασπορά στις τιμές της ορμής είναι μέγιστη
στο σημείο ισορροπίας, όπου έχουμε τη μέγιστη μέση ορμή, και ελάχιστη στη
μέγιστη απομάκρυνση, όπου η μέση ορμή είναι μηδέν. Στο αριστερό Σχήμα (1.8)
βλέπουμε στιγμιότυπα της ρ(x, t) για t ≈ T /8, T /4, T /2, 5T /8 και 3T /4.
0.8 0.8
t=0.785 t=0.785
t=1.571 t=1.571
0.7 0.7
t=3.142 t=3.142
t=3.927 t=3.927
0.6 t=4.712 0.6 t=4.712
0.5 0.5
ρ(x,t)
ρ(x,t)
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
-8 -6 -4 -2 0 2 4 6 8 -10 -5 0 5 10
x x
Σχήμα 1.8: Η πυκνότητα πιθανότητας ρ(x, t) για τον απλό αρμονικό ταλαντωτή, η οποία είναι
περιοδική συνάρτηση του χρόνου με περίοδο T = 2π/ω = 2π. Τα στιγμιότυπα είναι για t ≈
T /8, T /4, T /2, 5T /8 και 3T /4. (Αριστερά) Η ψ(x, 0) είναι η (1.25) με k0 = 3, σ = 0.5, x0 = 0.
(Δεξιά) Η ψ(x, 0) είναι η (1.25) με k0 = 6, σ = 0.5, x0 = 0.
Στη συνέχεια διπλασιάζουμε την k0 , την αρχική μέση ορμή της ψ(x, 0). Το
αποτέλεσμα το βλέπουμε στο δεξί Σχήμα (1.8). Το πλάτος της ταλάντωσης φαί-
νεται ότι έχει διπλασιαστεί, ενώ τα υπόλοιπα ποιοτικά χαρακτηριστικά της
ρ(x, t) παραμένουν τα ίδια.
Στη συνέχεια μπορούμε να√υπολογίσουμε τις αναμενόμενες τιμές ⟨x⟩, ⟨x2 ⟩
και την αβεβαιότητα ∆x = ⟨x2 ⟩ − ⟨x⟩2 , όπως στη σελίδα 17. Στο Σχήμα
1.9 βλέπουμε τα αποτελέσματα για τις συναρτήσεις ⟨x⟩(t) και ∆x(t). Η ⟨x⟩(t)
είναι ημιτονοειδής συνάρτηση του χρόνου με περίοδο T = 2π/ω = 2π και
πλάτους k0 . Η ∆x(t) είναι επίσης περιοδική συνάρτηση του χρόνου με περίοδο
T /2 = π. Παρατηρούμε ότι οι τιμές της είναι φραγμένες στο διάστημα [1/2, 1],
άρα το κυματοπακέτο δεν απλώνεται απεριόριστα, όπως στην περίπτωση του
ελεύθερου σωματιδίου.
Τα παραπάνω αποτελέσματα είναι αναμενόμενα. Οι συναρτήσεις ⟨x⟩(t) και
24 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
6 1
k0=3 k0=3
k0=6 0.95 k0=6
4
0.9
0.85
2
0.8
<x>(t)
Δx(t)
0 0.75
0.7
-2
0.65
0.6
-4
0.55
-6 0.5
0 1 2 3 4 5 6 0 1 2 3 4 5 6
t t
Σχήμα 1.9: (Αριστερά) Η αναμενόμενη τιμή ⟨x⟩(t) για τον απλό αρμονικό ταλαντωτή με αρχική
κατάσταση την (1.25) με σ = 0.5, x0 = 0. Η αναμενόμενη τιμή ⟨x⟩(t) είναι περιοδική με περίοδο
T = 2π/ω = 2π και με πλάτος k0 , και είναι ίδια με την απομάκρυνση x(t) του κλασικού
αρμονικού ταλαντωτή. (Δεξιά) Η αβεβαιότητα στη θέση ∆x(t) για τις δύο περιπτώσεις του
αριστερού σχήματος, η οποία είναι ανεξάρτητη του k0 .
∆x(t) υπολογίζονται στο Παράρτημα 1.8.3, όπου δείχνεται ότι είναι¹¹:
⟨x⟩(t) = k0 sin t
( ) ( )
1 1 1 1
2
(∆x(t)) = σ − 2 cos 2t +
2 2
σ + 2 . (1.54)
2 4σ 2 4σ
Η σύγκριση των τιμών που δείχνονται στο Σχήμα (1.9) με τις παραπάνω σχέσεις,
αφήνεται ως άσκηση στον αναγνώστη.
1.7 Μετρήσεις
Στις προηγούμενες ενότητες υπολογίσαμε αναμενόμενες τιμές φυσικών πο-
σοτήτων που είναι συνάρτηση της θέσης, όπως την ⟨x⟩∑ και την ⟨x2 ⟩. Το ολο-
κλήρωμα (1.11) προσεγγίστηκε από το άθροισμα ⟨x ⟩ ≈ i xai ρni ∆x και είδαμε
a
πως η προσέγγιση αυτή μας έδωσε αρκετά καλά αποτελέσματα. Στην ενότητα
αυτή θα υπολογίσουμε αναμενόμενες τιμές φυσικών ποσοτήτων που είναι συ-
ναρτήσεις και της ορμής, όπως λ.χ. τις ⟨p⟩, ⟨p2 ⟩ και ⟨H⟩. Επειδή p = −i∂/∂x,
υπολογισμός αυτός χρειάζεται εκτιμητές παραγώγων της κυματοσυνάρτησης
ψ(x, t).
Θα ακολουθήσουμε την προσέγγιση του Κεφαλαίου 10 και ο αναγνώστης θα
πρέπει να μελετήσει την ενότητα 10.4. Οι αλλαγές που χρειάζεται να κάνουμε
στο πρόγραμμα που βρίσκεται στο αρχείο observables.f90 είναι λίγες: Κυ-
ρίως πρέπει να λάβουμε υπόψη πως οι κυματοσυναρτήσεις που θα ολοκληρώ-
σουμε είναι μιγαδικές συναρτήσεις. Τον υπολογισμό των αναμενόμενων τιμών
θα τον απομονώσουμε σε μία subroutine
¹¹Δείτε τις σχέσεις (1.96), (1.113).
1.7. ΜΕΤΡΗΣΕΙΣ 25
subroutine observables ( psi , V , Nx , xmin , dx , avX , avP , avX2 ,&
avP2 , avE , DX2 , DP2 )
complex ( 8 ) , d i m e n s i o n ( Nx ) : : psi
...
στην οποία ο χρήστης της πρέπει να δώσει ένα μιγαδικό array psi(Nx), όπου
για δεδομένο tn , ψ(xi , tn ) → psi(i). Στην έξοδο παίρνει τις αναμενόμενες τι-
μές ⟨x⟩ → avX, ⟨p⟩ → avP, ⟨x2 ⟩ → avX2, ⟨p2 ⟩ → avP2, ⟨H⟩ → avE, (∆X)2 →
DX2, (∆P )2 → DP2. Παρατηρήστε πως στο πρόγραμμα αυτό, πρέπει να λά-
βουμε υπόψη ότι η κυματοσυνάρτηση παίρνει μιγαδικές τιμές, και ψ ∗ (xi , tn ) =
DCONJG(psi(i)). Επίσης παρατηρήστε πως, επειδή οι αναμενόμενες τιμές ερ-
μιτιανών τελεστών είναι πραγματικές, οι παραπάνω μεταβλητές ορίζονται να
είναι πραγματικές και το πρόγραμμα επιστρέφει το πραγματικό μέρος του ολο-
κληρώματος (1.11). Το πλήρες πρόγραμμα παρατίθεται παρακάτω και μπορεί
να βρεθεί στο αρχείο observables.f90 στο συνοδευτικό λογισμικό.
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
subroutine observables ( psi , V , Nx , xmin , dx , &
avX , avP , avX2 , avP2 , avE , DX2 , DP2 )
i m p l i c i t none
integer : : Nx
real (8) : : dx , xmin
complex ( 8 ) , d i m e n s i o n ( Nx ) : : psi
real ( 8 ) , d i m e n s i o n ( Nx ) : : V
complex ( 8 ) , a l l o c a t a b l e : : obs ( : ) , x ( : )
complex ( 8 ) : : integral
real (8) : : avX , avX2 , avP , avP2 , avV
real (8) : : avE , DX2 , DP2 , norm
integer :: i
complex ( 8 ) , p a r a m e t e r : : IMU = ( 0 . 0 D0 , 1 . 0 D0 )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
ALLOCATE ( obs ( Nx ) , x ( Nx ) )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! Define values of x :
do i = 1 , Nx
x ( i ) = xmin + ( i −1) * dx
end do
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! Norm o f p s i :
obs = psi * DCONJG ( psi )
c a l l integrate ( obs , dx , Nx , integral )
norm = DBLE ( integral )
psi = psi / norm
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! <x>
obs = psi * x * DCONJG ( psi )
26 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
c a l l integrate ( obs , dx , Nx , integral )
avX = DBLE ( integral )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! <p>
obs ( 1 ) = DCONJG ( psi ( 1 ) ) * ( psi ( 2 ) − psi ( 1 ) ) / dx
do i = 2 , Nx − 1
obs ( i ) = DCONJG ( psi ( i ) ) * ( psi ( i + 1 ) − psi ( i −1) ) / ( 2 . 0 D0 * dx )
end do
obs ( Nx ) = DCONJG ( psi ( Nx ) ) * ( psi ( Nx ) − psi ( Nx −1) ) / dx
c a l l integrate ( obs , dx , Nx , integral )
avP = DIMAG ( integral ) ! m u l t i p l y by − i
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! <x ^2 >
obs = psi * x * x * DCONJG ( psi )
c a l l integrate ( obs , dx , Nx , integral )
avX2= DBLE ( integral )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! <p^2 >
obs ( 1 ) = DCONJG ( psi ( 1 ) ) * ( 2 . 0 D0 * psi ( 1 ) −5.0 D0 * psi ( 2 )&
+ 4 . 0 D0 * psi ( 3 )−psi ( 4 ) ) / ( dx * dx )
do i = 2 , Nx−1
obs ( i ) = DCONJG ( psi ( i ) ) &
* ( psi ( i + 1 ) −2.0 D0 * psi ( i ) +psi ( i −1) ) / ( dx * dx )
enddo
obs ( Nx ) = DCONJG ( psi ( Nx ) ) * ( 2 . 0 D0 * psi ( Nx ) −5.0 D0 * psi ( Nx −1)&
+ 4 . 0 D0 * psi ( Nx −2)−psi ( Nx −3) ) / ( dx * dx )
c a l l integrate ( obs , dx , Nx , integral )
avP2=−DBLE ( integral )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! <V ( x ) >
obs = psi * V * DCONJG ( psi )
c a l l integrate ( obs , dx , Nx , integral )
avV = DBLE ( integral )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! <E>
avE = 0 . 5 D0 * avP2 + avV
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! U n c e r t a i n t i e s : DX2 , DP2 :
DX2 = avX2 − avX * avX
DP2 = avP2 − avP * avP
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
DEALLOCATE ( obs , x )
end s u b r o u t i n e observables
! =========================================================
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
subroutine integrate ( psi , dx , N , integral )
i m p l i c i t none
integer :: N
complex ( 8 ) : : psi ( N ) , integral , i n t
1.7. ΜΕΤΡΗΣΕΙΣ 27
real ( 8 ) : : dx
integer : : Nx , i
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! Nx must be odd :
Nx = N ; i f (MOD( N , 2 ) == 0 ) Nx = N−1
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! zeroth order point :
i = 1
i n t = psi ( i )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! odd o r d e r p o i n t s ( i =k +1 i s even ) :
do i = 2 , Nx −1 ,2
i n t = i n t + 4 . 0 D0 * psi ( i )
enddo
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! even o r d e r p o i n t s :
do i = 3 , Nx −2 ,2
i n t = i n t + 2 . 0 D0 * psi ( i )
enddo
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! l a s t point :
i = Nx
i n t = i n t + psi ( i )
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! measure n o r m a l i z a t i o n :
i n t = i n t * dx / 3 . 0 D0
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! I f N i s even , we add t h e l a s t p o i n t :
i f (MOD( N , 2 ) == 0 ) i n t = i n t + psi ( N ) * dx
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
!−−−−− f i n a l r e s u l t :
integral = i n t
end s u b r o u t i n e integrate
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
Στο κυρίως πρόγραμμα πρέπει να κάνουμε μικρές μεταβολές. Οι πιο σημα-
ντικές αφορούν τον ορισμό της κυματοσυνάρτησης. Αντί για τον ορισμό (1.37),
θα χρησιμοποιήσουμε τον ακριβέστερο ορισμό που δίνεται από τη σχέση (1.116)
i
ψin = Rin + (Iin + Iin−1 ) . (1.55)
2
Αυτό μπορεί να γίνει με την εντολή
psi=DCMPLX ( repsi ( new , : ) , 0 . 5 D0 * ( impsi ( new , : ) +impsi ( old , : ) ) )
όπου η σταθερά new δίνει τις συναρτήσεις στο νέο βήμα και η old στο προη-
γούμενο. Οι αλλαγές που πρέπει να γίνουν φαίνονται παρακάτω. Το πρόγραμμα
28 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
μπορεί να βρεθεί στο αρχείο tdse_fd_obs.f90 στο συνοδευτικό λογισμικό.
complex ( 8 ) , d i m e n s i o n ( P ) : : psi
real (8) : : avX , avP , avX2 , avP2 , avE , DX2 , DP2
...
open ( u n i t =f_obs , f i l e = ’ o b s ’ )
...
!−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
! B e g i n l o o p on t i m e e v o l u t i o n :
do it = 1 , Nt
....
i f (MOD( it , tskip ) == 0 ) t h e n ! make a measurement
...
psi=DCMPLX ( repsi ( new , : ) , 0 . 5 D0 * ( impsi ( new , : ) +impsi ( old , : ) ) )
c a l l observables ( psi , V , Nx , xmin , dx , &
avX , avP , avX2 , avP2 , avE , DX2 , DP2 )
w r i t e ( f_obs , 1 0 1 ) it , t , &
avX , avP , avX2 , avP2 , avE , DX2 , DP2
end i f
....
end do
Το πρόγραμμα μπορεί να μεταγλωτιστεί με την εντολή
> gfortran tdse_fd_obs . f90 observables . f90 −o t
και τρέχει ακριβώς όπως και το πρόγραμμα που έχουμε χρησιμοποιήσει μέχρι
τώρα. Τα αποτελέσματα βρίσκονται στο αρχείο obs και μπορείτε να τα δείτε με
το gnuplot:
gnuplot > plot ” obs ” using 2:3 w i t h lines t i t l e ” <x> ”
gnuplot > plot ” obs ” using 2:4 w i t h lines t i t l e ” <p> ”
gnuplot > plot ” obs ” using 2:5 w i t h lines t i t l e ” <x ^2 > ”
gnuplot > plot ” obs ” using 2:6 w i t h lines t i t l e ” <p^2 > ”
gnuplot > plot ” obs ” using 2:7 w i t h lines t i t l e ” <H> ”
gnuplot > plot ” obs ” using 2:8 w i t h lines t i t l e ” Dx^2 ”
gnuplot > plot ” obs ” using 2:9 w i t h lines t i t l e ” Dp^2 ”
gnuplot > plot ” obs ” using 2:( s q r t ( $8 * $9 ) ) w l t ” Dx Dp ”
Ας θεωρήσουμε πρώτα το σωμάτιο της ενότητας 1.4, και ειδικότερα αυτό
που μελετήθηκε στο Σχήμα 1.4. Το κυματοπακέτο που περιγράφει την κατά-
στασή του κινείται πρακτικά ελεύθερο για t ≲ 8, ενώ μετά αλληλεπιδρά με τα
τοιχώματα του πηγαδιού στις θέσεις x = ±xmax .
Για το ελεύθερο σωμάτιο, το οποίο βρίσκεται αρχικά στην κατάσταση (1.25),
1.7. ΜΕΤΡΗΣΕΙΣ 29
15 3
<x> <p>
2.5
10
2
5
1.5
0
1
-5
0.5
-10
0
-15 -0.5
-20 -1
-25 -1.5
0 5 10 15 20 25 30 0 5 10 15 20 25 30
t t
Δx Δp
16 3
14
12 2.5
10
8 2
6
1.5
4
2
1
0
0 5 10 15 20 25 30 0 5 10 15 20 25 30
t t
Σχήμα 1.10: Οι αναμενόμενες τιμές ⟨x⟩(t), ⟨p⟩(t) και οι αβεβαιότητες ∆x(t)∆p(t) για το ελεύ-
θερο σωμάτιο με x0 = −27, k0 = 3, σ = 0.5. Οι πράσινες καμπύλες υπολογίζονται από τις
σχέσεις (1.56), οι οποίες αποκλίνουν από τα αριθμητικά αποτελέσματα για t ≳ 8. Αυτό γίνεται
λόγω των ανακλάσεων του κυματοπακέτου στα τοιχώματα x = ±30.
η χρονική εξάρτηση των παρακάτω φυσικών ποσοτήτων υπολογίζεται¹² στο
Παράρτημα 1.8.3:
⟨x⟩(t) = k0 t + x0 ,
⟨p⟩(t) = k0 ,
( )
1 2 1 1
⟨H⟩(t) = ⟨p ⟩(t) = 2
k0 + 2 ,
2 2 4σ
√
1 2
∆x(t) = t + σ2 ,
4σ 2
1
∆p(t) = . (1.56)
2σ
Ειδικό ενδιαφέρον παρουσιάζει η ⟨H⟩(t), επειδή είναι διατηρούμενη ποσότητα
και μπορεί να χρησιμοποιηθεί ώς το μέγεθος που ελέγχει την ακρίβεια του αριθ-
μητικού υπολογισμού. Στο αριστερό Σχήμα 1.11 φαίνεται πως, μέσα στα πλαί-
σια της αριθμητικής ακρίβειας, η ⟨H⟩(t) παραμένει σταθερή. Στο Σχήμα 1.10
συγκρίνουμε τις φυσικές ποσότητες των (1.56) με τις αριθμητικά υπολογιζόμε-
νες. Παρατηρούμε πως οι τιμές ταυτίζονται μέχρι τα φαινόμενα της ανάκλασης
από τα τοιχώματα να γίνουν παρατηρήσιμα για t ≈ 8.
¹²Δείτε τις σχέσεις (1.88),(1.91),(1.93).
30 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
4.9979 60
<H> ΔX ΔP
4.99785
50
4.9978
40
4.99775
4.9977 30
4.99765
20
4.9976
10
4.99755
4.9975 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30
t t
Σχήμα 1.11: Η αναμενόμενη τιμή ⟨H⟩(t) και το γινόμενο των αβεβαιοτήτων ∆x(t)∆p(t) για
το ελεύθερο σωμάτιο με x0 = −27, k0 = 3, σ = 0.5. Η ⟨H⟩(t) = ⟨H⟩(0) = 5 διατηρείται
καθόλη τη διάρκεια της κίνησης του σωματιδίου. Η μικρή απόκλιση οφείλεται στο αριθμητικό
σφάλμα των υπολογισμών στο observables.f90, ενώ οι μικρές διακυμάνσεις στα αριθμητικά
σφάλματα της μεθόδου ολοκλήρωσης. Η αύξηση του ∆x(t)∆p(t) αντανακλά το γεγονός ότι
το κυματοπακέτο απλώνεται σε ολόκληρο το απειρόβαθο πηγάδι δυναμικού με τοιχώματα στα
x = ±30.
Στη συνέχεια θα μελετήσουμε τον αρμονικό ταλαντωτή, όπως κάναμε και
στην ενότητα 1.6, με x0 = 0, σ = 0.5, k0 = 3, ∆x = 0.02, a = 0.45, xmax 7.
Η χρονική εξάρτηση των παρακάτω φυσικών ποσοτήτων, υπολογίζεται¹³
στο Παράρτημα 1.8.3:
⟨x⟩(t) = k0 sin t ,
⟨p⟩(t) = k0 cos t ,
( )
1 2 1 1
⟨H⟩(t) = k0 + 2
σ + 2 ,
2 2 4σ
√ ( ) ( )
1 1 1
∆x(t) = √ + σ − 2 cos 2t + σ 2 +
2 ,
2 4σ 4σ 2
√ ( ) ( )
1 1 1
∆p(t) = √ − σ + 2 cos 2t + σ 2 +
2 . (1.57)
2 4σ 4σ 2
Στο Σχήμα 1.12 βλέπουμε τη χρονική εξάρτηση των ⟨x⟩(t), ⟨p⟩(t) και των αβε-
βαιοτήτων ∆x(t)∆p(t), ενώ στο Σχήμα 1.13 τη διατήρηση της ⟨H⟩(t) και το
γινόμενο των αβεβαιοτήτων ∆x(t)∆p(t). Παρατηρούμε πως το τελευταίο τα-
λαντώνεται με περίοδο T /4 μεταξύ μιας μέγιστης και μιας ελάχιστης τιμής. Η
συμφωνία των αριθμητικών υπολογισμών με τις σχέσεις (1.57) είναι εξαιρετική.
¹³Δείτε τις σχέσεις (1.96),(1.110),(1.113).
1.7. ΜΕΤΡΗΣΕΙΣ 31
3 3
<x> <p>
2 2
1 1
0 0
-1 -1
-2 -2
-3 -3
0 1 2 3 4 5 6 0 1 2 3 4 5 6
t t
1 1
Δx Δp
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
0 1 2 3 4 5 6 0 1 2 3 4 5 6
t t
Σχήμα 1.12: Οι αναμενόμενες τιμές ⟨x⟩(t), ⟨p⟩(t) και οι αβεβαιότητες ∆x(t)∆p(t) για τον αρ-
μονικό ταλαντωτή με x0 = 0, k0 = 3, σ = 0.5. Η θεωρητική πρόβλεψη που δίνεται από τη
σχέση (1.57) ταυτίζεται με τις καμπύλες που δίνουν οι αριθμητικοί υπολογισμοί, εκτός από την
περίπτωση του ∆p(t) όπου η πράσινη καμπύλη είναι η θεωριτική πρόβλεψη και η μωβ οι αριθ-
μητικοί υπολογισμοί.
5.122705 0.64
<H> ΔX ΔP
5.122704
0.62
5.122703
0.6
5.122702
0.58
5.122701
5.122700 0.56
5.122699
0.54
5.122698
0.52
5.122697
0.5
5.122696
5.122695 0.48
0 1 2 3 4 5 6 0 1 2 3 4 5 6
t t
Σχήμα 1.13: Η αναμενόμενη τιμή ⟨H⟩(t) και το γινόμενο των αβεβαιοτήτων ∆x(t)∆p(t) για για
τον αρμονικό ταλαντωτή με x0 = 0, k0 = 3, σ = 0.5. Η ⟨H⟩(t) = ⟨H⟩(0) = 5 διατηρείται
καθόλη τη διάρκεια της κίνησης του σωματιδίου. Η μικρή απόκλιση οφείλεται στο αριθμητικό
σφάλμα των υπολογισμών στο observables.f90, ενώ οι μικρές διακυμάνσεις στα αριθμητικά
σφάλματα της μεθόδου ολοκλήρωσης. Το γινόμενο ∆x(t)∆p(t) ταλαντώνεται και η τιμή του
μένει φραγμένη ανάμεσα στη μέγιστη και ελάχιστη τιμή που παίρνει στο χρόνο.
32 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
1.8 Παράρτημα
1.8.1 Απόδειξη Σχέσης (1.40)
Αναπτύσσουμε τη συνάρτηση R(x, t) κατά Taylor γύρω από το t
( ) ( )2
∆t ∆t ∂R(x, t) 1 ∆t ∂ 2 R(x, t)
R x, t + = R(x, t) + + + O(∆t3 )
2 2 ∂t 2! 2 ∂t2
( ) ( )2
∆t ∆t ∂R(x, t) 1 ∆t ∂ 2 R(x, t)
R x, t − = R(x, t) − + − O(∆t3 ) ,
2 2 ∂t 2! 2 ∂t2
(1.58)
και αφαιρούμε κατά μέλη:
( ) ( )
∆t ∆t ∂R(x, t)
R x, t + − R x, t − = ∆t + O(∆t3 ) . (1.59)
2 2 ∂t
Αντικαθιστούμε στο δεξί μέλος την δεύτερη εξίσωση από τις (1.33) και παίρ-
νουμε:
( ) ( )
∆t ∆t
R x, t + = R x, t − +
2 2
( )
1 ∂ 2 I(x, t)
+ ∆t − + V (x)I(x, t) + O(∆t3 ) . (1.60)
2 ∂x2
Στη συνέχεια αναπτύσσουμε τη συνάρτηση I(x, t) κατά Taylor γύρω από το x
∂I(x, t) 1 ∂ 2 I(x, t)
I (x + ∆x, t) = I(x, t) + ∆x + ∆x2 + O(∆x3 )
∂x 2! ∂x2
∂I(x, t) 1 ∂ 2 I(x, t)
I (x − ∆x, t) = I(x, t) − ∆x + ∆x2 − O(∆x3 ) , (1.61)
∂x 2! ∂x2
και προσθέτουμε κατά μέλη:
∂ 2 I(x, t)
I (x + ∆x, t) + I (x − ∆x, t) = 2I(x, t) + ∆x2 + O(∆x4 ) , (1.62)
∂x2
οπότε παίρνουμε
∂ 2 I(x, t) 1
2
= 2
(I (x + ∆x, t) − 2I(x, t) + I (x − ∆x, t)) + O(∆x2 ) . (1.63)
∂x ∆x
Αντικαθιστώντας στην (1.60) παίρνουμε
( )
∆t
R x, t + =
2
( )
∆t ∆t
= R x, t − − {I (x + ∆x, t) − 2I(x, t) + I (x − ∆x, t)}
2 2∆x2
+ ∆tV (x)I(x, t) + O(∆t3 ) + O(∆x2 ∆t) . (1.64)
1.8. ΠΑΡΑΡΤΗΜΑ 33
Επιλέγουμε xi = x, tn = t − ∆t/2,(οπότε παίρνουμε ) xi+1 =
( x + ∆x, ) xi−1 =
x − ∆x, tn+1/2 = t, tn = t + ∆t/2, R x, t + 2 = Ri , R x, t − 2 = Rin ,
∆t n+1 ∆t
I(x, t) = Iin , I (x + ∆x, t) = Ii+1
n
, I (x − ∆x, t) = Ii−1
n
, V (x) = V (xi ) ≡ Vi και
καταλήγουμε στη σχέση
( n )
Rin+1 = Rin − a Ii+1 − 2Iin + Ii−1
n
+ ∆tVi Iin + O(∆t3 ) + O(∆x2 ∆t) , (1.65)
όπου a ≡ ∆t/(2∆x2 ). Αυτή δεν είναι άλλη από την πρώτη εξίσωση στην (1.40),
αν θεωρήσουμε τα σφάλματα αμελητέα. Με τον ίδιο ακριβώς τρόπο αποδεικνύ-
εται και η δεύτερη.
Σύμφωνα με την (1.65), σε κάθε βήμα προστίθεται σφάλμα τάξης O(∆t3 ) και
O(∆x2 ∆t). Μετά από Nt βήματα, τα σφάλματα θα είναι
tmax
∼ Nt × O(∆t3 ) ∼ × O(∆t3 ) ∼ O(∆t2 ) , (1.66)
∆t
από τον πρώτο όρο, και
tmax
∼ Nt × O(∆x2 ∆t) ∼ × O(∆x2 ∆t) ∼ O(∆x2 ) , (1.67)
∆t
από τον δεύτερο. Αν κρατήσουμε το ∆x σταθερό, θεωρώντας το O(∆x2 ) αμε-
λητέο, τότε η μέθοδος δίνει σφάλμα τάξης ∆t2 , οπότε θεωρείται δεύτερης τάξης
ως προς ∆t.
Στην πράξη, όμως, θέλουμε να κρατήσουμε το ∆t όσο γίνεται μεγαλύτερο
και να μειώσουμε το συνολικό σφάλμα. Για ∆x αρκετά μικρό, (2/∆x)2 ≫ V+
και V+ + (2/∆x)2 ≈ 2/∆x)2 , οπότε η σχέση ευστάθειας (1.46), γίνεται
∆t 1
a= ≤ . (1.68)
2(∆x)2 2
Μειώνοντας το ∆t μαζί με το ∆x, κρατώντας το a σταθερό και κοντά στο 1,
έχουμε ∆x2 ∼ ∆t/(2a), οπότε το σφάλμα από τον όρο (1.67) γίνεται O(∆x2 )
∼ O(∆t) και επικρατεί από αυτό του όρου (1.66) και το συνολικό σφάλμα από
την (1.40) γίνεται O(∆t).
1.8.2 Απόδειξη Σχέσης (1.48)
Η πυκνότητα πιθανότητας εύρεσης του σωματιδίου στη θέση x δίνεται από
τη σχέση
ρ(x, t) = ψ(x, t)∗ ψ(x, t) = R(x, t)2 + I(x, t)2 . (1.69)
Η σχέση αυτή δεν μπορεί να χρησιμοποιηθεί με τον δυναμικό κανόνα (1.40),
γιατί τα R(x, t) και I(x, t) ορίζονται σε διαφορετικές χρονικές στιγμές. Οπότε
34 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
για τις χρονικές στιγμές t = n∆t ορίζουμε την ποσότητα:
∆t ∆t
)I(x, t −
ρ̂(x, t) = R(x, t)2 + I(x, t + )
2 2
∆t ∆t 2
ρ̂(x, t + ) = R(x, t)R(x, t + ∆t) + I(x, t + ) . (1.70)
2 2
Χρησιμοποιώντας τις σχέσεις (1.58) για τη συνάρτηση I(x, t ± ∆t
2
), παίρνουμε
∆t ∆t
I(x, t + )I(x, t − )=
( 2 2 )( )
∆t ∂I(x, t) ∆t ∂I(x, t)
= I(x, t) + + O(∆t ) 2
I(x, t) − + O(∆t ) =
2
2 ∂t 2 ∂t
( )
∆t ∂I(x, t) ∂I(x, t)
2
= I(x, t) + −I(x, t) + I(x, t) + O(∆t2 )) =
2 ∂t ∂t
= I(x, t)2 + O(∆t2 ) , (1.71)
άρα
ρ̂(x, t) = R(x, t)2 + I(x, t)2 + O(∆t2 ) = ρ(x, t) + O(∆t2 ) . (1.72)
Πρέπει να τονιστεί ότι το σφάλμα O(∆t2 ) της παραπάνω σχέσης δεν συσσω-
ρεύεται με τον χρόνο, αφού σε κάθε βήμα ο υπολογισμός γίνεται με τις νέες τιμές
των R(x, t)2 και I(x, t)2 . Στο σφάλμα του υπολογισμού της ρ(x, t) πρέπει κα-
νείς να προσθέσει τα σφάλματα των R(x, t)2 και I(x, t)2 . Με παρόμοιο τρόπο
μπορεί κανείς να δείξει ότι
∆t ∆t 2 ∆t 2
ρ̂(x, t + ) = R(x, t + ) + I(x, t + ) + O(∆t2 ) . (1.73)
2 2 2
Ορίζουμε
ρni = ρ̂(xi , tn )
∆t ∆t
= R(xi , tn )2 + I(xi , tn + )I(xi , tn − )
2 2
= (Rin )2 + Iin Iin−1
n+ 12 ∆t
ρi = ρ̂(xi , tn + )
2
∆t 2
= R(xi , tn )R(xi , tn + ∆t) + I(xi , tn + )
2
= Rin Rin+1 + (Iin )2 . (1.74)
1.8. ΠΑΡΑΡΤΗΜΑ 35
1
Η συνολική πιθανότητα P n και P n+ 2 που δίνεται από τις εξισώσεις (1.48)
στο σχήμα Visscher διατηρείται. Πράγματι, λαμβάνοντας υπόψη τις συνοριακές
συνθήκες
R1n = RN n
x
= I1n = INn x = 0 , (1.75)
έχουμε ότι
∑
N x −1 ∑
N x −1
[ n 2 ]
n
P = ρni = (Ri ) + Iin Iin−1
i=2 i=2
1
∑
N x −1
n+ 21
∑
N x −1
[ ]
P n+ 2 = ρi = (Iin )2 + Rin Rin+1 , (1.76)
i=2 i=2
οπότε
n+ 12
ρi − ρni = Rin (Rin+1 − Rin ) + Iin (Iin − Iin−1 )
[ ]
= Rin −a(Ii+1n
− 2Iin + Ii−1n
) + ∆tVi Rin Iin
[ ]
+ Iin +a(Ri+1n
− 2Rin + Ri−1 n
) − ∆tVi Iin Rin
= −aRin Ii+1
n
− aRin Ii−1
n
+ aIin Ri+1
n
+ aIin Ri−1
n
. (1.77)
Στη δεύτερη σειρά χρησιμοποιήσαμε τις εξισώσεις του σχήματος Visscher (1.40).
Αθροίζοντας ως προς i παίρνουμε
∑
Nx −1 ( ) ∑
Nx −1 ∑
Nx −1
n+ 1
ρi 2 − ρni = −a Rin Ii+1
n
−a Rin Ii−1
n
i=2 i=2 i=2
∑
N x −1 ∑
Nx −1
+a Iin Ri+1
n
+a Iin Ri−1
n
i=2 i=2
∑
Nx ∑
Nx −1
= −a n
Ri−1 Iin − a Rin Ii−1
n
i=3 i=2
∑Nx ∑
N x −1
n
+a Ii−1 Rin + a n
Iin Ri−1
i=3 i=2
= −aRNx −1 INx − aR2 I1
n n n n
n
+ aINn x −1 RN x
+ aI2n R1n
= 0, (1.78)
όπου στην τελευταία γραμμή χρησιμοποιήσαμε τις συνοριακές συνθήκες (1.75).
Άρα
1 1 1
P n+ 2 − P n = 0 ⇒ P n+ 2 = P n = P n− 2 = P n−1 = . . . = P 1 , (1.79)
36 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
και οι ποσότητες αυτές διατηρούνται ακριβώς.
1.8.3 Αναμενόμενες Τιμές
Για να υπολογίσουμε τις αναμενόμενες τιμές συναρτήσει του χρόνου (1.50),
(1.54), (1.56), (1.57), μπορούμε να χρησιμοποιήσουμε το θεώρημα του Ehrenfest:
d
i ⟨A(x, p)⟩(t) = ⟨[A(x, p), H]⟩ , (1.80)
dt
όπου A(x, p) είναι ένας ερμιτιανός τελεστής που δεν έχει ρητή εξάρτηση από
τον χρόνο, H είναι ο τελεστής της Χαμιλτονιανής και [A(x, p), H] = A(x, p)H−
HA(x, p) είναι ο μεταθέτης δύο τελεστών. Παρατηρήστε πως η παραπάνω εξί-
σωση είναι φορμαλιστικά όμοια με αυτή που προκύπτει από την εξίσωση κίνη-
σης της κλασικής μηχανικής
d
A(x, p) = {A(x, p), H}
dt
για την κλασική ποσότητα A(x, p), με την αντικατάσταση (1/i)[·, ·] → {·, ·},
όπου {A(x, p), H} είναι η αγγύλη Poisson¹⁴ των A(x, p), H.
Η απόδειξη του θεωρήματος είναι εύκολη. Έστω ότι το σωμάτιο βρίσκεται
στην κατάσταση ψ(x, t), οπότε:
∫
⟨A(x, p)⟩(t) = ψ(x, t)∗ A(x, p)ψ(x, t)dx . (1.81)
Παρατηρήστε πως η χρονική εξάρτηση της ⟨A(x, p)⟩(t) προκύπτει από τη χρο-
νική εξάρτηση της ψ(x, t), η οποία δίνεται από την εξίσωση του Schrödinger
∂ψ(x, t)
i = Hψ(x, t) . (1.82)
∂t
¹⁴{A, H} = (∂A/∂x)(∂H/∂p) − (∂H/∂x)(∂A/∂p)
1.8. ΠΑΡΑΡΤΗΜΑ 37
Τότε
d
⟨A(x, p)⟩(t) =
dt ∫ ∫
∂ψ(x, t)∗ ∂ψ(x, t)
= A(x, p)ψ(x, t)dx + ψ(x, t)∗ A(x, p) dx
∂t ∂t
∫ ∫
1 1
= ( Hψ(x, t))∗ A(x, p)ψ(x, t)dx + ψ(x, t)∗ A(x, p)( Hψ(x, t))dx
i i
∫ ∫
1 1
= (− )ψ(x, t)∗ H † A(x, p)ψ(x, t)dx + ψ(x, t)∗ A(x, p)H( )ψ(x, t)dx
i i
∫ ∫
1 1
=− ψ(x, t)∗ HA(x, p)ψ(x, t)dx + ψ(x, t)∗ A(x, p)Hψ(x, t)dx
i i
∫
1
= ψ(x, t)∗ (A(x, p)H − HA(x, p))ψ(x, t)dx
i
∫
1
= ψ(x, t)∗ [A(x, p), H]ψ(x, t)dx
i
1
= ⟨[A(x, p), H]⟩ ,
i
οπότε παίρνουμε την (1.81). Στην 3η γραμμή χρησιμοποιήσαμε τον ορισμό του
Ερμιτιανού συζυγούς¹⁵ H † της Χαμιλτονιανής και στην 4η το γεγονός ότι η Χα-
μιλτονιανή είναι ένας Ερμιτιανός τελεστής, δηλαδή H = H † .
Στους παρακάτω υπολογισμούς θα χρειαστούμε τη βασική σχέση τελεστών
ορμής-θέσης
[x, p] = i , (1.83)
και τις αλγεβρικές σχέσεις μεταξύ οποιωνδήποτε τελεστών
[A, B] = −[B, A] ,
[λA + µB, C] = λ[A, C] + µ[B, C] ,
[AB, C] = A[B, C] + [A, C]B , (1.84)
οι οποίες αποδεικνύονται απλά με ανάπτυγμα και των δύο πλευρών των εξισώ-
σεων.
Επίσης θα χρειαστούμε την αρχική τιμή των αναμενόμενων τιμών των παρα-
κάτω τελεστών, όταν το σωματίδιο βρίσκεται στην αρχική κατάσταση ψ(x, 0)
που δίνεται από την (1.25). Αυτές δίνονται από τις (1.27), (1.28), τις οποίες πα-
†
¹⁵Ο
∫ Ερμιτιανός ∫ ∗ † A του
τελεστής τελεστή
∫ A ορίζεται
∫ από τη σχέση ⟨ψ| A |ϕ⟩ = (⟨ϕ| A† |ψ⟩)∗
⇔ ψ Aϕ dx = ( ϕ A ψ dx) = ϕ(A ψ) dx = (A ψ)∗ ϕ dx. Προφανώς (A† )† = A.
∗ ∗ † ∗ †
38 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
ραθέτουμε παρακάτω:
⟨x⟩(0) = x0 , ⟨p⟩(0) = k0 ,
⟨x2 ⟩(0) = x20 + σ 2 , ⟨p2 ⟩(0) = k02 + 4σ1 2 ,
(1.85)
⟨xp⟩(0) = k0 x0 + 2i , ⟨px⟩(0) = k0 x0 − 2i ,
1
∆x(0) = σ , ∆p(0) = 2σ ,
οι οποίες αποδεικνύονται εύκολα κάνοντας τα ολοκληρώματα (1.11) με την
ψ(x, 0).
Έστω ένα ελεύθερο σωμάτιο με Χαμιλτονιανή
1
H = p2 , (1.86)
2
οπότε παίρνουμε
p2 p p p p
[x, H] = [x, ] = [x, p] + [x, p] = i + i = ip ,
2 2 2 2 2
2
p p p
[p, H] = [p, ] = [p, p] + [p, p] = 0 ,
2 2 2
[x2 , H] = x[x, H] + [x, H]x = i(xp + px) ,
[p2 , H] = p[p, H] + [p, H]p = 0 + 0 = 0 ,
[xp, H] = x[p, H] + [x, H]p = x0 + (ip)p = ip2 ,
[px, H] = [xp − i, H] = [xp, H] = ip2 . (1.87)
Εφαρμόζουμε την εξίσωση (1.80) για
d
i ⟨p⟩(t) = ⟨[p, H]⟩ = 0 ⇒ ⟨p⟩(t) = ⟨p⟩(0) = k0 , (1.88)
dt
d 1
i ⟨p2 ⟩(t) = ⟨[p2 , H]⟩ = 0 ⇒ ⟨p2 ⟩(t) = ⟨p2 ⟩(0) = k02 + 2 , (1.89)
dt 4σ
1
(∆p(t)) = (∆p(0)) = ⟨p ⟩(0) − (⟨p⟩(0)) = 2 .
2 2 2 2
(1.90)
4σ
και για
d
i ⟨x⟩(t) = ⟨[x, H]⟩ = i⟨p⟩(t) = i⟨p⟩(0) = ik0 ⇒
dt
⟨x⟩(t) = k0 t + ⟨x⟩(0) = k0 t + x0 , (1.91)
1.8. ΠΑΡΑΡΤΗΜΑ 39
d 2
i ⟨x ⟩(t) = ⟨[x2 , H]⟩ = i⟨xp + px⟩ ,
dt ( )
d 1
i ⟨xp⟩(t) = ⟨[xp, H]⟩ = i⟨p ⟩ = i k0 + 2 ⇒
2 2
dt 4σ
( ) ( )
1 1 i
⟨xp⟩(t) = k0 + 2 t + ⟨xp⟩(0) = k0 + 2 t + k0 x0 +
2 2
4σ 4σ 2
( ) ( )
1 1 i
⟨px⟩(t) = k0 + 2 t + ⟨px⟩(0) = k0 + 2 t + k0 x0 −
2 2
4σ 4σ 2
Αντικαθιστώντας τις δύο τελευταίες σχέσεις στην πρώτη, παίρνουμε
( )
d 2 1
⟨x ⟩(t) = ⟨xp⟩ + ⟨px⟩ = 2 k0 + 2 t + 2k0 x0 ⇒
2
dt 4σ
( )
1
⟨x2 ⟩(t) = k02 + 2 t2 + 2k0 x0 t + ⟨x2 ⟩(0)
4σ
( )
1
2
= k0 + 2 t2 + 2k0 x0 t + x20 + σ 2 . (1.92)
4σ
Η αβεβαιότητα ∆x υπολογίζεται από τις σχέσεις
(∆x(t))2 = ⟨x2 ⟩(t) − (⟨x⟩(t))2
{( ) }
1
= k0 + 2 t + 2k0 x0 t + x0 + σ − (k0 t + x0 )2
2 2 2 2
4σ
1 2
= 2 t + σ2 . (1.93)
4σ
Στη συνέχεια θεωρούμε τον απλό αρμονικό ταλαντωτή
1 1
H = p 2 + x2 . (1.94)
2 2
Για απλότητα, θα θεωρήσουμε το αρχικό κυματοπακέτο να έχει το κέντρο του
στην αρχή των αξόνων, δηλ. x0 = 0. Από τις σχέσεις
1
[x, H] = [x, p2 ] = ip ,
2
1 2 1 1 1 1
[p, H] = [p, x ] = x[p, x] + [p, x]x = x(−i) + (−i)x = −ix ,
2 2 2 2 2
παίρνουμε
d d2 d
i ⟨x⟩(t) = ⟨[x, H]⟩ = +i⟨p⟩ ⇒ 2
⟨x⟩(t) = + ⟨p⟩ , (1.95)
dt dt dt
2
d d d
i ⟨p⟩(t) = ⟨[p, H]⟩ = −i⟨x⟩ ⇒ 2
⟨p⟩(t) = − ⟨x⟩ ,
dt dt dt
40 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
οπότε
d2 d2
⟨x⟩(t) = −⟨x⟩ ⇒ ⟨x⟩(t) + ⟨x⟩ = 0 ,
dt2 dt2
d2 d2
⟨p⟩(t) = −⟨p⟩ ⇒ ⟨p⟩(t) + ⟨p⟩ = 0 .
dt2 dt2
Οι παραπάνω εξισώσεις δίνουν αρμονική ταλάντωση των μεγεθών και, λαμβά-
νοντας υπόψη τις αρχικές συνθήκες ⟨x⟩(0) = 0, ⟨p⟩(0) = k0 και την (1.95),
παίρνουμε
⟨x⟩(t) = k0 sin t ,
⟨p⟩(t) = k0 cos t . (1.96)
Για να υπολογίσουμε τις αναμενόμενες τιμές των x2 , p2 , H, είναι βολικό να
θεωρήσουμε τους τελεστές δημιουργίας και καταστροφής, a† και a αντίστοιχα:
1 1
a = √ (x + ip) , a† = √ (x − ip) ⇔ (1.97)
2 2
1 ( ) i ( )
x = √ a† + a , p = √ a† − a . (1.98)
2 2
από όπου προκύπτει η θεμελιώδης σχέση:
[ ]
† x + ip x − ip 1 1
[a, a ] = √ , √ = (−i[x, p] + i[p, x]) = (−ii + i(−i))
2 2 2 2
= 1. (1.99)
Ισχύουν οι σχέσεις:
( )2 ( )2
1 i † 1 1 †
H= √ (a − a) + √ (a + a)
2 2 2 2
1 ( )
= aa† + a† a
2
1
= a† a +
2
† 1
= aa − , (1.100)
2
1( † † )
2
x = +a a + aa† + a† a + aa , (1.101)
2
1( † † )
2
p = −a a + aa† + a† a − aa , (1.102)
2
1.8. ΠΑΡΑΡΤΗΜΑ 41
καθώς και οι
[a, a† ] = 1 , (1.103)
1
[a, H] = [a, a† a + ] = [a, a† ]a = a , (1.104)
2
1
[a† , H] = [a† , a† a + ] = a† [a† , a] = a† (−1) = −a† , (1.105)
2
1
[a2 , H] = [a2 , a† a + ] = [a, a† a]a + a[a, a† a] = 2a2 , (1.106)
2
1
[(a† )2 , H] = [(a† )2 , a† a + ] = [a† , a† a]a† + a† [a, a† a] = −2(a† )2 . (1.107)
2
Οπότε παίρνουμε
d i
i ⟨a⟩(t) = ⟨[a, H]⟩ = ⟨a⟩ ⇒ ⟨a⟩(t) = ⟨a⟩(0) e−it = √ k0 e−it ,
dt 2
i
⟨a† ⟩(t) = (⟨a⟩(t))∗ = − √ k0 eit ,
2
√ √
αφού ⟨a⟩(0) = (⟨x⟩(0) + i⟨p⟩(0))/ 2 = (0 + ik0 )/ 2.
d 2
i ⟨a ⟩(t) = ⟨[a2 , H]⟩ = 2⟨a2 ⟩ ⇒ ⟨a2 ⟩(t) = ⟨a2 ⟩(0) e−i2t ,
dt
⟨(a† )2 ⟩(t) = (⟨a2 ⟩(t))∗ = (⟨a2 ⟩(0))∗ ei2t .
Αλλά
⟨( )2 ⟩
x + ip
⟨a ⟩(0) =
2
√
2
1 1 i
= ⟨x2 ⟩(0) − ⟨p2 ⟩(0) + (⟨xp⟩ + ⟨px⟩)
2 ( 2 ) 2 ( )
1 2 1 1 i i i
= σ − 2
k0 + 2 + −
2 2 4σ 2 2 2
1 1 1
= σ 2 − k02 − 2 .
2 2 8σ
Οπότε τελικά:
( )
1 2 1 2 1
⟨a ⟩(t) =
2
σ − k − e−i2t , (1.108)
2 2 0 8σ 2
( )
† 2 1 2 1 2 1
⟨(a ) ⟩(t) = σ − k − e+i2t . (1.109)
2 2 0 8σ 2
42 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
Παρομοίως:
d † 1
i ⟨a a⟩(t) = ⟨[a† a, H]⟩ = ⟨H − , H⟩ = 0 ⇒ ⟨a† a⟩(t) = ⟨a† a⟩(0) .
dt 2
Αλλά
⟨ ⟩
† 1 2 1 2 1
⟨a a⟩(0) = x + p − (0) ,
2 2 2
( )
1 2 1 1 1
= σ + k0 + 2 − ,
2
2 2 4σ 2
οπότε λαμβάνοντας υπόψη ότι aa† = a† a + 1 και H = a† a + 1/2, παίρνουμε
( )
† 1 2 1 1 1
⟨a a⟩(t) = σ + k0 + 2 − ,
2
2 2 4σ 2
( )
1 1 1 1
⟨aa† ⟩(t) = σ 2 + k02 + 2 + ,
2 2 4σ 2
( )
1 1 1
⟨H⟩(t) = σ 2 + k02 + 2 . (1.110)
2 2 4σ
Επίσης
1( † † )
⟨x2 ⟩(t) = ⟨a a ⟩(t) + ⟨a† a⟩(t) + ⟨aa† ⟩(t) + ⟨aa⟩(t)
2 {( ) }
1 1 2 1 2 1 −i2t 1
= σ − k0 − 2 (e +i2t
+ e ) + σ + k0 + 2
2 2
2 2 2 8σ 4σ
( )
1 2 1 2 1 1 1 1
= σ − k0 − 2 cos 2t + σ 2 + k02 + 2 . (1.111)
2 2 8σ 2 2 8σ
Παρομοίως
1( )
⟨p2 ⟩(t) = −⟨a† a† ⟩(t) + ⟨a† a⟩(t) + ⟨aa† ⟩(t) − ⟨aa⟩(t)
2( )
1 2 1 2 1 1 1 1
=− σ − k0 − 2 cos 2t + σ 2 + k02 + 2 . (1.112)
2 2 8σ 2 2 8σ
Τέλος
( ) ( )
1 1 1 1
(∆x(t)) = ⟨x ⟩(t) − (⟨x⟩(t)) = +
2 2 2
σ − 2 cos 2t +
2 2
σ + 2 ,
2 4σ 2 4σ
( ) ( )
1 1 1 1
(∆p(t)) = ⟨p ⟩(t) − (⟨p⟩(t)) = −
2 2 2
σ + 2 cos 2t +
2 2
σ + 2 .
2 4σ 2 4σ
(1.113)
1.8. ΠΑΡΑΡΤΗΜΑ 43
Αφήνουμε σαν άσκηση να αποδείξετε ότι
( )
1 1 i
⟨xp⟩(t) = − σ 2 − k02 − sin 2t + ,
2 4σ 2 2
( )
1 1 i
⟨px⟩(t) = − σ 2 − k02 − sin 2t − , (1.114)
2 4σ 2 2
άρα ( )
1
⟨xp⟩(t) + ⟨px⟩(t) = − σ − k0 − 2 sin 2t
2 2
(1.115)
4σ
44 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
1.9 Ασκήσεις
1.1 Να αποδείξετε ότι η κυματοσυνάρτηση (1.25) είναι σωστά κανονικοποιη-
μένη και στη συνέχεια να αποδείξετε τις σχέσεις (1.27) και (1.28).
1.2 Να αποδείξετε την εξίσωση (1.73)
1.3 Η σχέση (1.35) υπολογίζει την κυματοσυνάρτηση με σφάλμα O(∆t). Ένας
πιο ακριβής ορισμός θα ήταν:
i
ψ(xi , tn ) = R(xi , tn ) + {I(xi , tn+ 1 ) + I(xi , tn− 1 )}
2 2 2
1
ψ(xi , tn+ 1 ) = {R(xi , tn+1 ) + R(xi , tn )} + iI(xi , tn+ 1 ) . (1.116)
2 2 2
Να αποδείξετε ότι το σφάλμα μειώνεται σε O(∆t2 ) και να υλοποιήσετε
τον ορισμό της ψ(xi , tn ) στα προγράμματά σας (αφορά μόνο τα προγράμ-
ματα που υπολογίζουν αναμενόμενες τιμές).
1.4 Στο πρόγραμμα tdse_fd.f90, θέσαμε την παράμετρο tinypsi = 1.0D-20.
Για να καταλάβετε το ρόλο της, να θέσετε tinypsi = 1.0D-300 και να
τρέξετε το πρόγραμμα για V (x) = 0, x0 = −7, σ = 0.5, k0 = 3, ∆x =
0.02, a = 0.45, xmax = 20, tmax = 0.3. Να κάνετε τη γραφική παρά-
σταση της πυκνότητας πιθανότητας ρ(x) για όλες τις τιμές του χρόνου t.
Τι παρατηρείτε; Γιατί; Να θέσετε πάλι tinypsi = 1.0D-20 και να επα-
ναλάβετε. Τι παρατηρείτε;
1.5 Το ελεύθερο σωμάτιο, το οποίο αρχικά βρίσκεται στην κατάσταση ψ(x, 0) =
(2/π)1/4 exp (−x2 + ik0 x), τη χρονική στιγμή t βρίσκεται στην κατάσταση
( ) { ( )2 }
(2/π)1/4 1 2 1 ik0
ψ(x, t) = √ exp − k0 exp − x− ,
1 + 2it 4 1 + 2it 2
(1.117)
με
{ }
(2/π)1/2 2(x − k0 t)2
ρ(x, t) = |ψ(x, t)| = √
2
exp − . (1.118)
1 + 4t2 1 + 4t2
Να αποδείξετε ότι η (1.117) είναι λύση της εξίσωσης Schrödinger και στη
συνέχεια να συγκρίνετε τα αποτελέσματα του προγράμματος tdse_fd.f90
με τις παραπάνω σχέσεις, επιλέγοντας k0 = 3, xmax = 40.
1.9. ΑΣΚΗΣΕΙΣ 45
1.6 Μελετήστε το δυναμικό της εξίσωσης (1.51) για V0 = 1, 2, 3, 4, 6, 8, 10,
12, 16, b = 0.5. Να θέσετε k0 = 3, σ = 0.5, ∆x =∫ 00.02, xmax = 30, tmax =
15. Να υπολογίσετε τις πιθανότητες pR (t) = −xmax ρ(x, t)dx, pT (t) =
∫ +xmax
0
ρ(x, t)dx και από αυτές τις pR,∞ = pR (+∞), pT,∞ = pT (+∞). Να
κάνετε τη γραφική παράσταση των pR,∞ , pT,∞ συναρτήσει του V0 και να
συζητήσετε τα αποτλέσματά σας.
1.7 Μελετήστε το δυναμικό της εξίσωσης (1.51) για V0 = 12, b = 0.01, 0.02,
0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.30, 0.40, 0.50. Να θέ-
σετε k0 = 3, σ = 0.5, ∆x = ∫ +xΝα
∫ 00.02, xmax = 30, tmax = 15. υπολογίσετε
max
τις πιθανότητες pR (t) = −xmax ρ(x, t)dx, pT (t) = 0 ρ(x, t)dx και
από αυτές τις pR,∞ = pR (+∞), pT,∞ = pT (+∞). Να κάνετε τη γραφική
παράσταση των pR,∞ , pT,∞ συναρτήσει του b και να συζητήσετε τα απο-
τλέσματά σας.
1.8 Να αναπαράξετε το Σχήμα 1.9 για k0 = 3, 6, 9, 12 και να συγκρίνετε τις
τιμές τους με το θεωρητικό αποτέλεσμα (1.54) κάνοντας τις γραφικές πα-
ραστάσεις (ti , ⟨x⟩(ti ) − ⟨x⟩theory (ti )), (ti , ∆x(ti ) − ∆xtheory (ti )), όπου οι
τιμές ⟨x⟩theory (t), ∆xtheory (t) θα υπολογίζονται από τις σχέσεις (1.54).
1.9 Να αναπαράξετε τα Σχήματα 1.12 και 1.13 παίρνοντας x0 = 0, k0 = 3,
σ = 0.5, ∆x = 0.02, a = 0.45, xmax = 12. Στη συνέχεια να κάνετε
τις γραφικές παραστάσεις των (ti , ⟨x⟩(ti ) − ⟨x⟩theory (ti )), (ti , ⟨x2 ⟩(ti ) −
⟨x2 ⟩theory (ti )), (ti , ⟨p⟩(ti )−⟨p⟩theory (ti )), (ti , ⟨p2 ⟩(ti )−⟨p2 ⟩theory (ti )), (ti , ∆x(ti )−
∆xtheory (ti )), (ti , ∆p(ti ) − ∆ptheory (ti )), όπου οι θεωρητικές τιμές θα υπο-
λογίζονται από τις σχέσεις (1.57). Να επαναλάβετε για k0 = 6, 12 επιλέ-
γοντας κατάλληλη τιμή για το xmax , έτσι ώστε η ψ(x, t) να είναι αμελητέα
στη γειτονιά των x = ±xmax
√
1.10 Coherent states: Παρατηρήστε πως όταν σ = 1/ 2 η σχέση (1.57) προ-
βλέπει πως η διασπορά ∆x(t) του κυματοπακέτου γίνεται ελάχιστη και
ανεξάρτητη του χρόνου. Το τελευταίο σημαίνει πως το σχήμα της ρ(x, t)
δεν παραμορφώνεται με τον χρόνο. Το σωμάτιο αυτό βρίσκεται σε μια κα-
τάσταση η οποία μπορεί να θεωρηθεί πως είναι η “κοντινότερη” σε αυτή
ενός κλασικού αρμονικού ταλαντωτή: Μία coherent state ακολουθεί την
κλασική τροχιά του σωματιδίου όσο γίνεται κοντύτερα.
(αʹ) Να επαναλάβετε
√ την προηγούμενη άσκηση για x0 = 0, k0 = 3,
σ = 1/ 2, ∆x = 0.02, a = 0.45, xmax = 12, tmax = 30. Πόσο είναι
το αριθμητικό σφάλμα στη ∆x(t);
(βʹ) Να επαναλάβετε το προηγούμενο ερώτημα για ∆x = 0.01
46 ΚΕΦΑΛΑΙΟ 1. ΚΒΑΝΤΙΚΗ ΔΥΝΑΜΙΚΗ
(γʹ) Σε κάθε μια από τις προηγούμενες περιπτώσεις, να κάνετε τη γρα-
φική παράσταση της ρ(x, t) για χρόνους tm ≈ mT /8.
(δʹ) Να επαναλάβετε τα παραπάνω για k0 = 6, xmax = 15.
(εʹ) Να διαβάσετε για τις coherent states στη Wikipedia
https://en.wikipedia.org/wiki/Coherent_states
και να εμβαθύνετε περισσότερο τη μελέτη σας πάνω σε αυτές. Να
επιβεβαιώσετε ότι η ψ(x, 0) που χρησιμοποιήσατε είναι ειδική περί-
πτωση της γενικότερης μορφής των κυματοσυναρτήσεων ψ (α) (x, t)
που αναφέρονται στην παραπάνω ιστοσελίδα και να αντιστοιχή-
σετε τις παραμέτρους στην ψ (α) (x, t) με αυτές που χρησιμοποιήθη-
καν στην ψ(x, 0) (εξίσωση (1.25)).
ΒΙΒΛΙΟΓΡΑΦΙΑ
[Κεφάλαιο 1]
[1] P.B. Visscher, “A fast explicit algorithm for the time–dependent Schrödinger
equation”, Computers in Physics 5 (1991) 596.
47
48 ΒΙΒΛΙΟΓΡΑΦΙΑ