Compressibilities and Phonon dispersions of cluster crystals

(1)

Compressibilities and Phonon dispersions of cluster crystals

Diplomarbeit Johannes H¨ aring

1. Gutachter (Betreuer): Herr Professor Dr. M. Fuchs 2. Gutachter: Herr Professor Dr. G. Maret

12.12.2013

(2)

(3)

Bibliography

[1] C. N. Likos et al., J. Chem. Phys. 126, 224502 (2007) [2] C. Walz and M. Fuchs, Phys. Rev. B 81, 134110 (2010)

[3] T. Neuhaus and C. N. Likos, J. Phys: Condens. Matter 23 (2011) 234112

[4] A. P. Gantapara, Soft Condensed Matter Physics Group, Debye Institute, Utrecht, work in progress

[5] B. M. Mladek, P. Charbonneau and Daan Frenkel, PRL 99, 235702 (2007) [6] C. Walz, G. Szamel and M. Fuchs, unpublished

[7] K. Zhang, Patrick Charbonneau and B. M Mladek, PRL 105, 245701 (2010) [8] N. G. Wilding and P. Sollich, EPL, 101 (2013) 10004

[9] Christof Walz, Dissertation, Ver¨offentlicht im Konstanzer Online-Publikations- System (KOPS) der Universit¨at Konstanz

[10] B. M. Mladek et al., PRL 96, 045701 (2006)

(5)

1 Einleitung

In dieser Arbeit sollen grundlegende Eigenschaften eines einfachen Modells eines Clus- terkristalls numerisch berechnet werden. Unter Clusterkristall ist wie beim normalen Kristall eine periodische Anordnung von Gitterpl¨atzen gemeint, die aber auch von mehr- eren Teilchen besetzt werden k¨onnen. Zur theoretischen Verwirklichung dieser Situation wird das GEM-4-Modell verwendet. GEM-4 bedeutet generalized exponential model mit dem Exponenten vier. Konkret bedeutet das ein isotropisches Paarpotential der Form

Φ(r) =˜ e⁻⁽^σ^r⁾⁴ (1.1)

mit dem Abstand r, der Reichweite des Potentials σ und der maximalen Energie . Gegenüber einem Hartkugelmodell ist das Potential durchnach oben beschränkt. Dies erlaubt mehrere Teilchen pro Gitterplatz. Es mag zuerst überaschen, dass bei einem rein repulsiven Potential überhaupt eine Clusterbildung stattfindet. Und tatsächlich tritt die Clusterbildung bei GEM-m-Potentalen mit m <3 nicht auf. Für die Untersuchung der Eigenschaften des Modells wird die mittlere Teilchendichte ρ aber einfach als gegeben angenommen. Das System kann sich also nur entscheiden, ob höhere Dichten über eine Verringerung der Gitterkonstante a oder eine Erhöhung der mittleren Besetzungszahl n_c erreicht werden. Ob ein System zu letzterem neigt, hängt vom Vorzeichen der Fouri- ertransformierten ab. Hat diese auch negative Anteile, ist Clusterbildung möglich. Die Rolle der negativen Anteile der Fouriertransformierten des GEM-4-Potentials wird später klar werden.

In [10] wurden mit Hilfe der Dichtefunktionaltheorie (DFT) die Phasen des Systems, die Werte der Gitterkonstanten und mittleren Besetzungszahl bestimmt. Dies soll nochmals nachvollzogen werden. Als Ansatz f¨ur die Teilchendichte in den festen Phasen wird auch hier verwendet:

n(~r) = nc(α π)³² X

R~

e^−α(~^r−^R)^~ ² (1.2)

mit der mittleren Besetzungszahl n_c, den Gittervektoren R~ und dem Quadrat der inversen Breite der Dichteverteilung um einen Gitterplatz α. Die Art des Gitters (FCC oder BCC) wird ¨uber die Summation ¨uber die Gittervektoren festgelegt. Durch den Ansatz wird aus dem Freien Energie-Funktional eine Funktion, die bei gegebener mit- tlerer Dichte und Temperatur nach (α,n_c) minimiert werden kann. Wie man sehen wird, ergibt dies neben der fluiden Phase sowohl eine Kristallisation im BCC- als auch im FCC-Gitter.

Mit den DFT-Optimalwerten von (α,n_c) ist die Dichteverteilung n(~r) vollständig bestimmt. Dies ermöglicht die Bestimmung weiterer Eigenschaften des Systems. Aus einer Arbeit von C. Walz und M. Fuchs ([6]) ergibt sich damit die Möglichkeit, über

(6)

eine Wellengleichung die Dispersionsrelationen bei verschiedenen mittleren Dichten und Temperaturen zu berechnen. Der verwendete Formalismus geht auf die Dissertation von C. Walz zurück ([9]) und verbindet phänomenologische Ausdrücke der verallgemein- erten Hydrodynamik des Kristalls mit mikroskopischen Ausdrücken. Die Besonderheit des Formalismus ist, dass der Schwerpunkt nicht auf idealen Kristallen liegt, sondern ebenso Punktdefekte berücksichtigt werden. Dies schlägt sich in folgendem Ansatz für die Teilchendichtevariationen δn(~r) nieder

δn(~r) =−ρ∇ ·~u(~r)−δc(~r) (1.3) mit der mittleren Teilchendichte ρ, dem Verschiebungsfeld ~u(~r) und der Defektdichte (in negativer Notation) c(~r). Das Verschiebungsfeld stellt die Verzerrungen der Gitter- struktur dar und ist beim idealen Kristall die einzige Quelle f¨ur die Dichtefluktuationen.

Für das GEM-4-Modell und seine Clusterbildung bietet sich somit folgende Verknüpfung an: Bei Dichteänderungen entspricht das Verschiebungsfeld einer Veränderung der Git- terkonstanten und die Defektfluktuationen einer Schwankung der Besetzungszahl. Da die Besetzungszahl des Clusterkristalls durchaus Werte über zehn annehmen kann und vor allem im Prinzip variabel ist, sollten sich Schwankungen der Besetzungszahl eher deutlich machen als bei normalen Kristallen ohne Clusterbildung. Mit der hier erwähnten Beset- zungszahl ist stets die mittlere Besetzungszahl pro Gitterplatz gemeint, die somit kontinuierliche Werte annehmen kann. In der Realität ist natürlich jeder Gitterplatz nur mit einer ganzahligen Annzahl an Teilchen besetzt. Damit die Mittelung über die Teilchen- zahl Sinn ergibt, muss deswegen eine Bewegung der Teilchen zwischen den Gitterplätzen möglich sein. Oder anders ausgedrückt: Die Schwankungen der Besetzungszahl müssen groß genug sein. Tatsächlich zeigen Simulationen des GEM-4-Modells bei Temperatur gegen null nur ganzzahlige Besetzungszahlen und somit isostruktuelle Phasenübergänge.

Der Fall tiefer Temperaturen und der ¨Ubergang ganzzahliger zu kontinuierlicher Beset- zungszahl soll deswegen ebenfalls ber¨ucksichtigt werden.

Zusätzlich zur Dispersionrelation, die für alle Wellenvektoren innerhalb der ersten Bril- louin -Zone berechnet werden kann, soll noch speziell der Limes kleiner Wellenvektoren, also großer Wellenlängen betrachtet werden. Dies führt zu Größen wie der isothermen Kompressibilität und der Schallgeschwindigkeit.

(7)

2 Theorie

2.1 Die Freie Energie

Die genaue Herleitung des Freien Energie-Funktionals ist in [1] zu finden. Die Freie Energie ist aufgeteilt in einen idealen Anteil F_id[n] und einem Wechselwirkungsanteil F_ex[n]. Dieser stellt eine funktionale Taylorentwichlung dar. Die Funktionalableitungen sind proportional zu den direkten Korrelationsfunktionen. Nach der zweiten Ordnung kann die Entwicklung abgebrochen werden (Ramankrishnan-Yussouff-Näherung). Als Ansatz für die Zweiteilchen-Direkte-Korrelationsfunktionc(r~₁, ~r₂) wird eine Mean-Field- Näherung (MFA) vorgenommen:

c(~r₁, ~r₂)≈ −βΦ(~˜ r₂−~r₁) =−βΦ(r)˜ (2.1) mit dem Potential des GEM-4-Modells ˜Φ und dem Abstandsbetrag r = |r2 −r1|. Ab jetzt soll der dimensionslose Ausdruck Φ(r) = ˜Φ(r) benutzt werden.

F[n] =F_id[n] +F_ex[n] (2.2) Fid[n] = 1

β Z

d³~r[n(~r) ln[n(~r)Λ³]−n(~r)] (2.3) F_ex[n] =

2 Z

d³~r₁n(~r₁) Z

d³~r₂Φ(~r₁−~r₂)n(~r₂) (2.4) (2.5) Als erstes soll die Freie Energie für n(~r) = ρ = const., also für das Fluid berechnet werden. ρ= ^N_V ist dann einfach die gesamte Teilchenzahl N durch das Gesamtvolumen V. FürF_id^fl gilt dann

F_id^fl =F_id[ρ] = 1 β

Z

d³~rρ[ln(ρΛ³)−1] (2.6)

= 1

β[ln(ρΛ³)−1]

Z

d³~rρ= N

β[ln(ρΛ³)−1] (2.7)

Das Doppelintegral von F_ex^fl kann durch die Substitution~r⁰ :=~r₁ −~r₂ in zwei einzelene Integrale zerlegt werden

F_ex^fl =F_ex[ρ] = ρ² 2

Z d³~r₁

Z

d³~r⁰Φ(~r⁰) = N ρ 2

Z

d³~r⁰Φ(~r⁰) (2.8) Durch eine weitere Substitution ~r := ^~^r

0

σ wird das Integral ¨uber Φ dimensionslos und somit numerisch berechenbar.

Z

d³~r⁰Φ(~r⁰) =σ³ Z

d³~re^−r⁴ =:σ³Φ₀ (2.9)

(8)

F_ex^fl und F_id^fl zusammen ergeben F^fl = N

β[ln(ρΛ³)−1] +N ρσ³

2 Φ₀ (2.10)

Um auf dimensionslose Größen zu kommen, wird der Ausdruck für die Freie Energie durch N geteilt. Da diese Größen fest geben sind, ändern sie an der Minimierung nichts. Mit den Defintionent^∗ := _β¹ = ^k^B^T, ρ^∗ :=ρσ³ und f := _N^F ergibt sich

f^fl=t^∗(lnρ^∗+ 3 lnΛ

σ −1) + ρ^∗

2 Φ₀ (2.11)

Für die festen Phasen wählen wir folgenden Ansatz für die Dichte n(~r) =n_c(α

π)³² X

R~

e^−α(~^r−^R)^~ ² (2.12)

• n_cist die durchschnittlicheAnzahl Teilchen pro Gitterplatz, also ein kontinuierlicher Wert. Da die Teilchenzahl in Wirklichkeit ein diskreter Wert ist, wird sich zeigen m¨ussen, wann dieser Ansatz gerechtfertigt ist.

• P

R~ ist die Summation ¨uber alle Gitterpl¨atzeR. Hier fließt ein, welcher Gittertyp~ angenommen wird.

• e^−α(~^r−^R)^~ ² ist ein Ansatz f¨ur die Dichteverteilung um einen Gitterplatz.

• α ist ein Maß f¨ur das inverse Quadrat der Breite der Dichtverteilung.

Die Normierungskonstante (^α_π)³² stellt schließlich sicher, dass gilt Z

d³~rn(~r) = n_cX

R~

Z

d³~r(α

π)³²e^−α(~^r−^R)^~ ² =n_cX

R~

1 =N (2.13)

In einem Kristall sind die Teilchen recht stark um ihren Gitterplatz lokalisiert, bzw tritt nach dem Lindemannkriterium bei nennenswerter Auslenkung (ca. 10% der Git- terkonstante) schon der Schmelzpunkt ein. Somit ist der ¨Uberlapp der einzelnen Gauss- verteilungen minimal und es kann folgende N¨aherung vorgenommen werden

n(~r) ln[n(~r)] =n_c(α π)³² X

R~

e^−α(~^r−^R)^~ ²ln[n_c(α π)³² X

R~

e^−α(~^r−^R)^~ ²] (2.14)

≈n_c(α π)³² X

R~

e^−α(~^r−^R)^~ ²[ln[n_c(α

π)³²]−α(~r−R)~ ²] (2.15) Damit ergibt sich f¨ur den idealen Anteil der Freien Energie des Kristalls

F_id^solid ≈ 1 β

X

R~

n_c Z

d³~r(α

π)³²e^−α(~^r−^R)^~ ²[ln[n_c(α

π)³²Λ³]−α(~r−R)~ ²−1] (2.16)

(9)

Mit P

R~n_c=N und der Substitution~r⁰ =~r−R~ erh¨alt man F_id^solid = N

β[(lnnc+3 2lnα

π + 3 ln Λ−1) Z

d³~r⁰(α

π)³²e^−α~^r⁰² −

− Z

d³~r⁰(α

π)³²α~r⁰²e^−α~^r⁰²] (2.17) Das erste Integral ist duch die Normierung eins und f¨ur das zweite gilt

Z

d³~r⁰(α

π)³²e^−α~^r⁰² = ( Z

dr⁰(α

π)¹²e^−αr⁰²)³ = 1 (2.18) Z

d³~r⁰(α

π)³²α~r⁰²e^−α~^r⁰² = Z

d³~r⁰(α

π)³²α(r_x⁰²+r⁰²_y +r_z⁰²)e^−α~^r⁰²

= 3 Z

d³~r⁰(α

π)³²αr⁰²_xe^−α~^r⁰² = 3 Z

dr⁰_x(α

π)³²r⁰²_xe^−αr⁰²^x Z

dr_y⁰dr⁰_ze^−α(r⁰²^y^+r⁰²^z⁾ = 3

2 (2.19) Im letzten Schritt wird bei vorderen Integral partiell integriert um auf den Wert ¹₂ zu kommen. Mit α^∗ =ασ² folgt

F_id^solid = N

β(lnn_c+ 3

2lnα^∗− 3

2lnπ+ 3 lnΛ σ −5

2) (2.20)

Die Berechnung von F_id^solid vereinfacht sich stark, wenn das zweite Integral als Faltung angesehen wird

f(~r₁) :=

Z

d³~r₂n(~r₂)Φ(~r₁−~r₂) (2.21) Mit der Identit¨at (diese folgt durch simples Ausrechnen mit der Fouriertransformation)

Z

d³~rg(~r)f(~r) =

Z d³~q

(2π)³g^∗_~_qf_~_q (2.22) und der Tatsache, dass die Fouriertransformation eines Faltungsintegrals dem Produkt der einzeln transformierten Funktionen entspricht, ergibt sich

F_ex^solid[n] = 2

Z d³~q

(2π)³n^∗_~_qn_~_qΦ_~_q (2.23) Dies ist erst einmal keine Vereinfachung des Problems. Durch die Fouriertransformation der einzelnen Funktionen kann dies aber erreicht werden.

n_~_q =n_c(α π)³² X

R~

Z

d³~re^−i~^q~^re^−α(~^r−^R)^~ ² (2.24)

~r→~r+R~i

= n_c(α π)³² X

R~

Z

d³~re^−i~^q(~^r+^R)^~ e^−α~^r² (2.25)

=nc(α π)³² X

R~

e^−i~^{q ~}^R Z

drxe^−iq^x^r^xe^−αr²^x Z

drye^−iq^y^r^ye^−αr²^y Z

drze^−iq^z^r^ze^−αr²^z (2.26)

=n_ce⁻^q

2x 4αe⁻

q2 y 4αe⁻^q

2z 4α

X

R~

e^−i~^{q ~}^R=n_ce⁻^q

2 4α

X

R~

e^−i~^{q ~}^R (2.27)

(10)

Hierfür wird die Tatsache genutzt, dass die Fouriertransformierte einer Gaußverteilung wieder eine Gaußverteilung ist. Die Summe über alle Gitterplätze kann als dreidimen- sionale Fourierreihe mit den Koeffizienten gleich eins angesehen werden

n_~_q=n_ce⁻^q

2 4α

X

G~

(2π)³ 1 Va

δ(~q−G)~ (2.28)

=ρe⁻^q

2

4α X

G~

(2π)³δ(~q−G)~ (2.29)

(2.30) V_aist das Volumen pro Gitterplatz und somit die inverse Dichte der Gitterpl¨atzeρ_l = _n^ρ

c. G~ sind die Vektoren des reziproken Gitters, die die Bedingung e^{−i ~}^{G ~}^R = 1 ∀R~ erf¨ullen.

Da Φ(r) rotationssymmetrisch ist, ist auch Φ_q rotationssymmetrisch und f¨ur den dritten Faktor von Gl. 2.23 ergibt sich

Φ_~_q = Φ_q = 4π Z

drre^−iqr Z ∞

r

dr⁰r⁰Φ(r⁰) (2.31) Setzt man nun beide Ausdr¨ucke 2.27 und 2.29 mit 2.31 in 2.23 ein, folgt

F_ex^solid[n] =ρn_cX

G~

X

R~

Z

d³~qe⁻^q

2

2αe^−i~^{q ~}^Rδ(~q−G)Φ~ _~_q (2.32)

=ρn_cX

G~

e⁻^G

2 2α

X

R~

e^{−i ~}^{G ~}^RΦ_G_~ (2.33)

= N ρ 2

X

G~

e⁻^G

2

2αΦ_G_~ (2.34)

Schließlich muss das Ergebnis noch in dimensionslosen Größen ausgedrückt werden. Sub- stituiert man in der Fouriertransformation Φ_q mit~r⁰ =~r/σ erhält man

Φ_~_q=σ³ Z

d³~r⁰e^−i~^qσ~^r

0

e^−r

04

=σ³Φ_~_qσ (2.35)

Zusammen mit K~ =Gσ~ folgt

F_ex^solid[n] = N ρ 2 σ³X

G~

e⁻^G

2σ2

2ασ2ΦGσ~ (2.36)

= N ρ^∗ 2

X

K~

e⁻^K

2

2α∗Φ_K_~ (2.37)

(2.38) Das ergibt f¨ur die gesamte (dimensionslose) Freie Energie pro Teilchen

f^solid = F_id^solid +F_ex^solid

N (2.39)

=t^∗(lnn_c+3

2lnα^∗− 3

2lnπ+ 3 lnΛ σ − 5

2) + ρ^∗ 2

X

K~

e⁻^K

2

2α∗Φ_K_~ (2.40)

(11)

Es gilt f¨ur dieK-Vektoren mit der Gitterkonstante a und den dimensionslosen Basisvek-~ toren des reziproken Gitters~g_i

K~ =Gσ~ = 2πσ

a(n1~g1+n2~g2+n3~g3) (2.41) Dien_igehen von−∞bis∞und jedem (n₁, n₂, n₃)-Tripel ist genau ein Punkt im reizipro- ken Gitter zugeordnet. Die ~g_i sind nur vom Gittertyp abh¨angig. Die Summe P

K~ geht also in eine Summe P

(n1,n2,n3) uber. F¨¨ ur die freie Energie relevant sind also nur die (dimensionslose) Gitterkonstante ^a_σ und die Wahl des Gittertyps. Sei z die Anzahl der Gitterpl¨atze pro koventioneller Einheitszelle mit dem Volumen a³ und N_a die Anzahl der konventionellen Einheitszellen. Dann gilt

a³ σ³ = V

N_a N_a

N zn_c

σ³ = zn_c

ρ^∗ (2.42)

Es gibt also bei gegebenen ρ^∗ und Gittertyp (mit zB z = 4 für fcc und z = 2 für bcc eine klare Zuordnung zwischen (dimensionsloser) Gitterkonstante und Besetzungzahln_c. Ist der Gittertyp und das (t^∗, ρ^∗)-Paar gegeben, lässt sich die Freie Energie also nach (a^∗, n_c) oder (a^∗,â_σ) minimieren. Welche Phase bei gegebenen (t^∗, ρ^∗) vorliegt, muss durch Vergleich der absoluten Werte der jeweils minimierten Freien Energie-Funktion getesten werden. Aus drei Gründen ist es vorteilhaft f^fl von f^solid abzuziehen

∆f =f^solid−f^fl =t^∗(lnn_c ρ^∗ +3

2lnα^∗ −3

2lnπ− 3 2) + ρ^∗

2 X

K6=0~

e⁻^K

2

2α∗ΦK~ (2.43) (2.44) Erstens ist je nach Vorzeichen von ∆f klar welche Phase vorliegt, zweitens fällt der un- bekannte Term ^Λ_σ raus. Da f^fl nicht von (a^∗, n_c) abhängt, ändert sich am Minimum der Funktion nichts. Drittens kann dann durch Teilen durch ρ^∗ das gesamte Phasenprob- lem enorm vereinfacht werden, da weder die Minimierung von ∆f noch das Vorzeichen dadurch beinflusst wird (Die Dichte ist immer größer Null)

f˜(n_c

ρ^∗, α^∗, t^∗

ρ^∗) = ∆f ρ^∗ = t^∗

ρ^∗(lnn_c ρ^∗ +3

2lnα^∗− 3

2lnπ− 3 2) + 1

2 X

K6=0~

e⁻^K

2

2α∗Φ_K_~ (2.45) Es gibt also nur noch eine thermodynamische Zustandsgröße _ρ^t^∗∗. Die Phasenübergänge im (t^∗, ρ^∗)-Diagramm liegen alle auf Strahlen vom Nullpunkt aus. Das alles aber unter der Bedingung, dass n_c kontinuierliche Werte annehmen kann.

2.2 Λ

_ab

(~ q)

Die Λ_ab(~q)-Matrix ist in [2] Gl. (26) durch die Wellengleichung definiert und wird ben¨otigt um die Dispersionsrelation zu berechnen

Λ_ab(~q) = λ_ab(~q)−iq_aµ_b(~q) +iq_bµ^∗_a(~q) +q_aq_bν(~q) (2.46)

(12)

Die drei Funktionen sind gegeben durch ν(~q) = 1

βV Z

d³~r₁ Z

d³~r₂n(~r₁)n(~r₂)[δ(~r₁−~r₂)

n(~r1) −c(~r₁−~r₂)]e^−i~^q(~^r¹^−~^r²⁾ (2.47) µ_a(~q) = 1

βV Z

d³~r₁ Z

d³~r₂n(~r₁)∇_an(~r₂)c(~r₁−~r₂)(1−e^−i~^q(~^r¹^−~^r²⁾) (2.48) λab(~q) = 1

βV Z

d³~r1

Z

d³~r2∇an(~r1)∇bn(~r2)c(~r1−~r2)(1−e^−i~^q(~^r¹^−~^r²⁾) (2.49) Nun sollen die Funktionen im Einzelnen f¨ur das gegebene GEM4-Problem berechnet werden. Da die Korrelationsfunktion c(~r) durch den MFA-Ansatz proportional zum Potiential ist, ist sie ebenfalls rotationssymmetrisch. Der erste Summand von ν(~q) ist einfach

1 βV

Z d³~r1

Z

d³~r2n(~r1)n(~r2)δ(~r₁−~r₂)

n(~r₁) e^−i~^q(~^r¹^−~^r²⁾ (2.50)

= 1 βV

Z d³~r₁

Z

d³~r₂n(~r₂)δ(~r₁−~r₂)e^−i~^q(~^r¹^−~^r²⁾ (2.51)

= 1 βV

Z

d³~r₁n(~r₁) = N βV = ρ

β (2.52)

Der zweite Summand weist starke Ähnlichkeiten mit dem Ausdruck für den Wechsel- wirkungsanteil der Freien Energie F_ex[n] (Gl. 2.4) auf. Der Lösungsweg ist deshalb analog. Etwas Vorsicht ist bei der ~q-Abhängigkeit der Funktionen geboten. Darum erfolgt die Intergration nach~q⁰. Für die Faltung wird diesmal definiert

f(~r₁) :=

Z

d³~r₂n(~r₂)c(~r₁−~r₂)e^−i~^q(~^r¹^−~^r²⁾ (2.53)

=−β Z

d³~r₂n(~r₂)Φ(~r₁−~r₂)e^−i~^q(~^r¹^−~^r²⁾ (2.54) Die Fouriertransformierte des~r=~r1−~r2-Teils ist

Φ_~_q+~_q⁰ = Z

d³~rΦ(~r)e^−i(~^q+~^q⁰^)~^r (2.55) Dadurch erh¨alt man

βν(~q) =ρ− 1 V

Z d³~r1

Z

d³~r2n(~r1)n(~r2)c(~r1−~r2)e^−i~^q(~^r¹^−~^r²⁾ (2.56)

=ρ− 1 V

Z

d³~rn(~r)f(~r) (2.57)

=ρ− 1 V

Z d³~q⁰

(2π)³n^∗_~_q0f_~_q⁰ (2.58)

=ρ+ β V

Z d³~q⁰

(2π)³n^∗_~_q0n_~_q⁰Φ_~_q⁰_+~_q (2.59) (2.60)

(13)

Das weitere Vorgehen erfolgt dann analog zu Gl.(2.23) ff. Es wird also (2.27) und (2.29) in (2.59) eingesetzt

βν(~q) =ρ+βρn_c V

X

R~

X

G~

Z

d³~q⁰e⁻^q

02

2αe^−i~^q⁰^R^~δ(~q⁰−G)Φ~ _~_q⁰_+~_q (2.61)

=ρ+βρ²X

G~

e⁻^G

2

2αΦG+~~ q (2.62)

Da t^∗ = _β¹ gilt, kommt man schließlich auf ν(~q)

=t^∗ρ+ρ²X

G~

e⁻^G

2

2αΦ_G+~_~ _q (2.63)

F¨ur die Herleitung von Gl.(2.48) und Gl.(2.49) soll die Korrelationsfunktion c(~r) gleich durch −βΦ(~r) ersetzt und die Gleichungen durch geteilt werden

µ_a(~q)

= 1

V Z

d³~r₁ Z

d³~r₂n(~r₁)∇_an(~r₂)Φ(~r₁−~r₂)(e^−iq(~^r¹^−~^r²⁾−1) (2.64) λ_ab(~q)

= 1

V Z

d³~r₁ Z

d³~r₂∇_an(~r₁)∇_bn(~r₂)Φ(~r₁−~r₂)(e^−iq(~^r¹^−~^r²⁾−1) (2.65) Um das gleiche Prinzip wie bei ν anwenden zu k¨onnen, braucht man die Fouriertrans- formierte der Ableitung von n(~r). Dies erfolgt ¨uber die partielle Integration

F T[∇αn(~r)]~q = Z

d³~re^−~^q~^r∇αn(~r) (2.66)

= 0− Z

d³~r(−i)q_αe^−~^q~^rn(~r) (2.67)

=iq_αn_~_q (2.68)

(2.69) Das bedeutet f¨ur µ_a(~q)

µ_a(~q)

= 1

V

Z d³~q⁰

(2π)³n^∗_~_q⁰(iq⁰_a)n~q⁰(Φ~q⁰+~q−Φ~q⁰) (2.70)

=iρnc

V X

R~

X

G~

Z

d³~q⁰q_a⁰e⁻^q

02

2αe^−i~^q⁰^R^~δ(~q⁰ −G)(Φ~ _~_q⁰_+~_q−Φ_~_q⁰) (2.71)

=iρ²X

G~

G_ae⁻^G

2

2α(Φ_~_q⁰_+~_q−Φ_~_q⁰) (2.72) (2.73)

(14)

Beiλ_ab(~q) muss beachtet werden, dass eine Fouriertransformierte komplex konjugiert ist λab(~q)

= 1

V

Z d³~q⁰

(2π)³(−iq_a⁰)n^∗_~_q⁰(iq_b⁰)n_~_q⁰(Φ_~_q⁰_+~_q−Φ_~_q⁰) (2.74)

= 1 V

Z d³~q⁰

(2π)³q_a⁰q⁰_bn^∗_~_q0n_~_q⁰(Φ_~_q⁰_+~_q−Φ_~_q⁰) (2.75)

=ρn_c V

X

R~

X

G~

Z

d³~q⁰q_a⁰q⁰_be⁻^q

02

2αe^−i~^q⁰^R^~δ(~q⁰−G)(Φ~ _~_q⁰_+~_q−Φ_~_q⁰) (2.76)

=ρ²X

G~

G_aG_be⁻^G

2

2α(c_G_~ −c_G+~_~ _q) (2.77) Nun l¨asst sich Λ_ab komplett angeben

Λ_ab

(~q) = q_aq_bρ (2.78) +ρ²X

G~

e⁻^G

2

2α[(q_aG_b+q_bG_a+G_aG_b)(Φ_~_q⁰_+~_q−Φ_~_q⁰) +q_aq_bΦG~⁰+~q]

Mit~k =~qσ,K~ =Gσ~ und Φ_q =σ³Φ_k und multipliziert mit ^σ_ρ² ergibt sich die dimensionslose Form

σ²

ρΛ_ab(~k) = k_ak_bt^∗ +ρ^∗X

K~

e⁻^K

2

2α∗[(kaKb+kbKa+KaKb)(Φ_|K+~ ~k|+ ΦK) +kakbΦ_|K+~ ~k|] (2.79)

2.3 Λ

_ab

(~ q → 0)

Im Limes großer Wellenl¨angen gilt nach [2]

ν(~q)≈ν(0) +O(q²) (2.80) µ_a(~q)≈iµ_abq_b+O(q²) (2.81) µab = 1

βV Z

d³~r1

Z

d³~r2n(~r1)∇an(~r2)c(~r1−~r2)(~r1−~r2)b (2.82) (2.83) Im Grunde kann µ_ab gleich behandelt werden wie µ_a(~q), also im Fourierraum mit Hilfe des Faltungsintegrals intergriert werden. Um die Anzahl der Integrale minimal zu halten muss aber erst die Rotationssymmetrie desr~₁−r~₂-Integranden wiederhergestellt werden.

Mitr~₁−r~₂ =~r ergibt sich Z

d³~re^−i~^q~^rc(~r)r_b = Z

d³~ri ∂

∂qb

e^−i~^q~^rc(~r) (2.84)

=i ∂

∂q_bc~q =i∂q

∂q_b

∂

∂qcq = qb

qic⁰_q (2.85)

(15)

mit q6= 0 und

ic⁰_q =i ∂

∂qc_q = 4π Z

drre^−iqr Z ∞

r

dr⁰r⁰c(r⁰) (2.86) Der Fall q = 0 spielt wie man sehen wird keine Rolle. Einfach betrachtet funktion- iert die Wiederherstellung der Rotationssymmetrie des Integranden über die Ersetzung r_b = ^q_q^br. Die umständliche Herleitung über die Ableitung der Fouriertransformation macht aber deutlich, dass es sich nicht um eine simple Ersetzung durch eine beliebige Funktion~q(~r) = q(r)^~^r_r handelt, die man wegen ihrer explizieten~r-Abhängigkeit natürlich nicht einfach vor das Integral ziehen darf. Durch den zusätzlichen Faktor r ist der zu fouriertransformierende Term antisymmetrisch. Dadurch sind lediglich die imaginären Sinusanteile der Fourierkoeffizienten ungleich Null.

ic⁰_q=−i4π Z

drrsin(−qr) Z ∞

r

dr⁰r⁰c(r⁰) (2.87)

=i4π Z

drrsin(qr) Z ∞

r

dr⁰r⁰c(r⁰) = −iβ4π Z

drrsin(qr) Z ∞

r

dr⁰r⁰Φ(r⁰) (2.88)

=:−iβΦ⁽¹⁾_q (2.89) Φ⁽¹⁾q ist also reellwertig und schon hinsichtlich der numerischen Rechnung definiert. Der Rest l¨auft analog zur Herleitung von µ_a(~q)

µ_ab = 1

V

Z d³~q

(2π)³n~q^∗(iqa)n~q

q_b

q(−iΦ⁽¹⁾_q ) (2.90)

=ρnc

V X

R~

X

G~

Z

d³~qq_ae⁻^q

2

2αe^−i~^{q ~}^Rδ(~q−G)~ qb

qΦ⁽¹⁾_q (2.91)

=ρ²X

G6=0~

e⁻^G

2

2αG_aG_j,b

G Φ⁽¹⁾_G (2.92)

Nun soll kurz der Fall q = 0 besprochen werden. Die linke Seite von Gl.(2.84) lautet dann durch die Antisymmetrie von c(r)r_b bez¨uglich der b-Komponente

Z

d³~rc(~r)r_b = 0 (2.93)

Der q →0-Limes vonλ_ab(~q)≈λ_abcdq_cq_d+O(q⁴) lautet mitr~₁−r~₂ =~r λ_abcd = 1

β2V Z

d³~r₁ Z

d³~r₂∇_an(~r₁)∇_bn(~r₂)c(~r)r_cr_d (2.94) F¨ur die Wiederherstellung der Symmetrie muss also auf zwei Ableitungen zur¨uckgegriffen werden. Wieder mit q 6= 0:

Z

d³~re^−i~^q~^rc(~r)r_cr_d=− Z

d³~r ∂

∂qc

∂

∂qd

e^−i~^q~^rc(~r) (2.95)

=− ∂

∂q_c

∂

∂q_dc_~_q =−∂q

∂q_c

∂q

∂q_d( ∂

∂q)²c_q =−qc

q qd

qc”_q (2.96)

= 4π Z

drr²e^−iqr Z ∞

r

dr⁰r⁰c(r⁰) (2.97)

(16)

Bei q = 0 und für c = d verschwindet das Integral diesmal nicht, bleibt aber wegen der e^−r⁴-Abhängigkeit endlich. Da im Ausdruck für λ_abcd mit G_a und G_b multipliziert wird und diese bei G~₀ = 0 null sind, kann der Summand trotzdem ignoriert werden.

Der Integrand von (2.97) ist wegen r² symmetrisch und deswegen bleibt nur der reelle Kosinusanteil ¨ubrig

4π Z

drr²e^−iqr Z ∞

r

dr⁰r⁰c(r⁰) = 4π Z

drr²cos(−qr) Z ∞

r

dr⁰r⁰c(r⁰) (2.98)

= 4π Z

drr²cos(qr) Z ∞

r

dr⁰r⁰c(r⁰) =: −βΦ⁽²⁾_q (2.99) Damit erhalten wir schließlich

λ_abcd

=−1 2ρ²X

G6=0~

e⁻^G

2

2αG_cG_dG_cG_d

G² Φ⁽²⁾_G (2.100)

(17)

3 Resultate

3.1 Fouriertransformierte des Potentials

Die genauen Eigenschaften des Systems werden letztendlich durch die direkte Korre- lationsfunktion c(r) bestimmt. Da diese durch den MFA-Ansatz proportional zum Potential (c(r) = −βΦ(r) = −βe⁻⁽^r^σ⁾⁴) ist und hier ausschließlich dessen Fourier- transformierte verwendet wird, lohnt sich ein Blick auf diese. In Abb. 3.1 sind die wesentlichen Eigenschaften schnell zu erkennen. Erstens geht der Wert schnell gegen null für große, wodurch in der numerischen Berechnung Summationen bis zum dritten oder viertern Nachbarn im reziproken Raum ausreichend sind, zweitens dominieren die Werte nahe null und drittens ist ein klares Minimum beiqσ = 5.57 zu erkennen. Da der Wechselwirkungsanteil des Freien Energie-Funktionals minimiert wird,wird das System versuchen möglichst starke Beiträge zur Freien Energie nahe des Minimums zu haben.

Zur Erinnerung:

F_ex^solid[n] = N ρ^∗ 2

X

K~

e⁻^K

2

2α∗Φ_K_~ (3.1)

Nach Gl. 2.45 wird bei gebebenem (t^∗, ρ^∗) durch α^∗ und _σâ bzw ⁿ_ρ^c minimiert. Betra- chtet man die obige Summe, wird man annehmen, dass der erste Nachbar möglichst im Minimum liegt, also K1 = G1σ = 2π^σ_ad ≈ 5.5 mit der dimensionslosen und nur vom Gittertyp abhängigen Konstante d, die dem Betrag eines dimensionslosen Basisvektor des reziproken Gitters entspricht, also dem Wert d_{f cc} = √

3 f¨ur ein FCC-Gitter und dbcc =√

2 f¨ur ein BCC-Gitter. Dies entspricht einem Wert von ^a_σ ≈ 1.954 (FCC) bzw.

a

σ ≈1.595 (BCC), und zwar ohne (direkte) Ab¨angigkeit von (t^∗, ρ^∗). Da sich ^a_σ direkt in

nc

ρ umrechnenen lässt, neigt das System also dazu höhere Dichten fast nur über steigende nc zu erreichen und die Gitterkonstante quasi konstant zu lassen.

(18)

-0.5 0 0.5 1 1.5 2 2.5 3 3.5 4

0 2 4 6 8 10 12

Φqσ

qσ

-0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25

4 5 6 7 8 9 10 11 12

Φqσ

qσ

Figure 3.1: Die Fouriertransformierte des GEM4-Potentials. Das untere Bild ist ein Ausschnitt des oberen um das relevante Minimum besser sichtbar zu machen

(19)

3.2 Verhalten der Freien Energie

In diesem Abschnitt soll das genaue Verhalten der Freien Energie im FCC-Gitter betrachtet werden. Dazu wird jeweils einer der drei Parameter (ⁿ_ρ∗^c, α^∗,_ρ^t^∗∗) fixiert und die anderen beiden variiert. Wie aus dem vorigen Abschnitt bekannt, sollte die Git- terkonstante weitgehend konstant sein. Abb. 3.2 zeigt den Verlauf der Freien Energie bei fixierter Gitterkonstante. Die durchgezogene Linie entspricht dem t^∗ = 0-Fall, und damit einer freien Energie ohne idealem Anteil. Bez¨uglich α^∗ gibt es kein Minimum, f wird mit zunehmendem α^∗ immer kleiner. Der Faktor e⁻^K

2

2α∗ geht gegen eins. Der (negative) Beitrag des ersten Nachbarn wird also maximiert. Die stärkeren positiven Beiträge weiter entfernt liegender Nachbarn gleichen das nicht mehr aus, da ihre Werte betragsmäßig zu klein sind. Bei _ρ^t^∗∗ > 0 ändert sich das Bild. Der ideale Anteil von Gl.

2.45 skaliert logarithmisch mitα^∗. Deshalb sinkt der Optimalwert vonα^∗ mit steigendem

t^∗

ρ^∗. Die Verteilung der Teilchen um einen Gitterplatz wird bei steigender Temperatur also immer breiter, was den Erwartungen an thermische Bewegungen entspricht. Abb.

3.3 zeigt die Freie Energie bei fixiertem α^∗ und bestätigt die unwesentliche Änderung der Gitterkonstante. In Abb. 3.4 ist der Verlauf für verschiedene inverse α^∗ zu sehen, also für verschiedene Gewichtungene⁻^K

2

2α∗. Bei hohen inversenα^∗, also mit zunehmender Nachbargeneration sinkender Gewichtungen, steigt die Gitterkonstante tendenziell. Eine steigende Gitterkonstante bedeutet aber im reziproken Raum ein kleineres qσ. Große qσ sorgen im Grunde dafür, dass weiter entferntere Nachbargenerationen einen Beitrag möglichst nahe null liefern. Sind diese aber durch sinkende Gewichtungen unbedeu- tend, wird der Beitrag der nahen Nachbarn wichtig. Obwohl steigene Gitterkonstanten (kleinere qσ) die Lage des nächsten Nachbarn aus dem Minimum des Potentials bewe- gen, scheint dies trotzdem günstiger zu sein. Der Gewinn durch sinkende Beiträge der

¨

ubernächsten Nachbarn muss diesen Nachteil aufwiegen. Zusammenfassend lässt sich also sagen, dass die Gitterkonstante im Großen und Ganzen unveränderlich ist. Da aber steigende Temperaturen einen deutlichen Einfluss auf α^∗ haben, das die Gewichtung der einzelnen Beiträge der Nachbargenerationen bestimmt, kommt so eine kleine Variation der Gitterkonstante und somit auch von ⁿ_ρ∗^c durch die ”Hintertür”. Der Verlauf der Kurven in der Nähe der Minima sieht allgemein unkompliziert aus und sollte deswegen numerisch gut zu behandeln sein.

(20)

-0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3

50 100 150 200 250 300 350 400

∆f/ρ*

ασ² 0

0.05 0.1 0.15

Figure 3.2: Die normierte Freie Energie ohne DFT-Optimierung f¨ur vier verschiedene _ρ^t^∗∗. Die Gitterkonstante ^a_σ wurde fixiert

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6

1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4

∆f/ρ*

a/σ 0 0.5 1.0 1.5

Figure 3.3: Die normierte Freie Energie ohne DFT-Optimierung f¨ur vier verschiedenet^∗ bei ρ^∗ = 10. Das Quadrat der inversen Breite der Gaußverteilungen um die Gitterpunkte ασ² wurde fixiert

.

(21)

-0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4

∆/ρ*

a/σ 1/α=0 1/α=0.1 1/α=0.2 1/α=0.3

Figure 3.4: Die normierte Freie Energie ohne DFT-Optimierung f¨ur vier verschiedene Breiten _ασ¹2 beit^∗ = 0

3.3 Minimierung der Freien Energie und die Phasen¨ uberg¨ ange des Systems

3.3.1 Der Bereich kontinuierlicher Besetzungszahlen (t

^∗

& 0.1)

Durch den vorigen Abschnitt ist die grundlegende Situation der Optimalwerte zur Min- imierung der Freien Energie bekannt, was gute Anfangswerte zur numerischen Min- imierung der Freien Energie liefert. Die Minimierung wird jeweils für das BCC und das FCC-Gitter mit Hilfe von ˜f(ⁿ_ρ∗^c, α^∗,_ρ^t^∗∗) ( siehe Gl. 2.45) vorgenommen. Ist ˜f kleiner null, liegt die fluide Phase vor. Für den Übergang BCC ↔ FCC werden die jeweils minimierten Funktionen verglichen, die in Abb. 3.6 oben zu sehen sind. Auffallend ist der ausgesprochen kleine Unterschied zwischen den beiden festen Phasen (genauer zu sehen in Abb. 3.6 unten), der um fast zwei Größenordnungen kleiner als der Unter- schied zum Fluid. Vergleiche zur Simulation ([5]) zeigen auch deswegen eine signifikante Abweichung im Wert von _ρ^t^∗∗ des BCC ↔ FCC- Übergangs. Die Unterschiede in der Freien Energie sind dann aber trotzdem noch klein. Der Übergang BCC↔Fluid stimmt dagegen gut mit dem oberen Bereich der aus der Simulation berechneten Phasenkoex- istenzen der BCC/Fluid-Phase überein. Die Phasenübergänge FCC ↔ BCC ↔ Fluid stehen aber nicht im Fokus der weiteren Arbeit, da dort die FCC-Phase im Vordergrund steht. Wichtig sind vor allem die (ⁿ_ρ∗^c, α^∗)-Optimalwerte für gegebenes (_ρ^t^∗∗). Unter der Bedingung kontinuierlicher Besetzungszahlen sind (ⁿ_ρ∗^c, α^∗) also Funktionen von (_ρ^t^∗∗). In

(22)

Abb. 3.5 und 3.7 sehen wir diese Funktionen für FCC und BCC. Die Werte für die Gitterkonstante liegen in der FCC-Phase ca. bei zwei und in der BCC-Phase ca bei 1.6, also nahe den erwarteten Werten nach Abschnitt 3.1. Da sie leicht drüber und vor allem mit steigendem (_ρ^t^∗∗) größer werden, bestätigt die Überlegungen aus Abschnitt 3.2. Der Schmelzpunkt tritt ca. beiα^∗ = 20 auf. Durch den Wert der dimensionslosen Gitterkon- stante _σâ lässt sich α^∗ =ασ² in αa² umrechnen. Dessen inverse Wurzel entspricht dann der maximalen Auslenkung der Teilchen um den jeweiligen Gitterplatz, es ergibt sich α^−0.5 ≈0.14a und entspricht damit dem Lindemannkriterium. Abb. 3.8 und 3.9 zeigen noch einzelne Werte der FCC-Phase und deren Fits. Die sehr gute Übereinstimmung mit den Fits in Abb. 3.9 ist der schwachen Änderung der Werte und der damit verbun- denen guten Möglichkeit zur Taylorentwichlung geschuldet. Bei α^∗ dagegen dominiert die Skalierung mit (_ρ^t^∗∗)⁻¹, was der Skalierung des Auslenkungsquadrats der Teilchen mit (_ρ^t^∗∗) entspricht. Sei <∆x∆x > das (mittlere) Auslenkungsquadrat. Dann gilt:

<∆x∆x >=at^∗

ρ^∗ +b(t^∗

ρ^∗)² = t^∗

ρ^∗(a+bt^∗

ρ^∗) (3.2)

Multipliziert mit n_c ergibt sich

n_c<∆x∆x >=t^∗n_c

ρ^∗(a+bt^∗

ρ^∗)) (3.3)

Da ⁿ_ρ∗^c quasi konstant und b klein ist entspricht dies weitgehend dem Gleichvertei- lungssatz.

Compressibilities and Phonon dispersions of cluster crystals