How crystals form: A theory of nucleation pathways

Teoria

Każda teoria nukleacji wymaga dwóch elementów: po pierwsze, sposobu obliczania własności klastrów, takich jak ich struktura i energie swobodne, a po drugie, dynamicznego opisu fluktuacji. Klasyczna teoria funkcjonałów gęstości (cDFT) (21, 22) od dawna uznawana jest za potencjalnie zdolną do dostarczenia pierwszego elementu. Od czasu pionierskiej pracy Oxtoby’ego i Evansa w latach 80-tych (23), cDFT była używana do określania struktury i energii klastrów krytycznych dla, po pierwsze, zarodkowania kropel cieczy z pary, a później dla krystalizacji. Wczesne obliczenia zawierały wiele upraszczających założeń, ale nowsze prace wykazały ilościową zgodność z symulacją (24). Zaletą cDFT w porównaniu z alternatywnymi technikami, takimi jak teoria kryształu z polem fazowym (PFC) (25) i modele interfejsu dyfuzyjnego jest to, że jest ona w zasadzie „ab initio”, wymagając jedynie potencjału międzyatomowego jako danych wejściowych i że daje ilościowo dokładny opis korelacji i struktury w skali molekularnej. cDFT jest podstawową teorią, z której inne mogą być rozumiane jako przybliżenia (22).

W używanej tutaj implementacji cDFT, lokalna gęstość liczbowa każdego gatunku chemicznego jest dyskretyzowana na sześciennej siatce obliczeniowej, z odstępem siatki znacznie mniejszym niż poszczególne molekuły. Równowagowe lokalne gęstości są wyznaczane poprzez minimalizację funkcji lokalnych gęstości, dając zarówno równowagowy rozkład molekuł, jak i energię swobodną układu. Teoria poprawnie opisuje cechy w skali molekularnej, takie jak upakowanie gęstych płynów w warstwy w pobliżu ściany. W szczególności, podczas gdy jednorodna ciecz ma jednorodną gęstość, ciało stałe jest z natury niejednorodne w skali molekularnej, ponieważ gęstość jest ostro szczytowa w miejscach sieci i przechodzi do bardzo niskich wartości pomiędzy nimi. Ostatnie postępy w cDFT rozszerzyły jej zastosowanie na wysoce niejednorodne układy, takie jak gęste krople cieczy i stałe skupiska w równowadze z tłem pary o niskiej gęstości (patrz Rys. 1).

Drugim elementem niezbędnym do opisania zarodkowania jest opis fluktuacji (26). Naturalną ramą do tego jest hydrodynamika fluktuacyjna (FH), która została zapoczątkowana przez Landaua i była intensywnie badana i rozwijana. FH jest obecnie szeroko używanym narzędziem, które zostało zastosowane w wielu dziedzinach, takich jak teoria sprzężenia trybów, przejście szkliste i zarodkowanie, a jej podstawy w bardziej fundamentalnej mechanice statystycznej zostały ustalone. Podstawowymi wielkościami wykorzystywanymi w teorii są przestrzennie zmienne lokalne gęstości poszczególnych gatunków, jak również pola prędkości i temperatury. Dla dużych cząsteczek, takich jak koloidy lub makromolekuły, w kąpieli z mniejszych cząsteczek (np. wody), można wyprowadzić przybliżony efektywny opis dużych cząsteczek, w którym wpływ mniejszych cząsteczek jest modelowany jako kombinacja tarcia i siły stochastycznej. Jeśli tłumienie kąpieli jest silne, to można to sprowadzić do pojedynczego równania opisującego pole gęstości dla dużych cząsteczek (28, 29) mającego postać∂∂tn^t(r)=D∇⋅n^t(r)∇δFδn^t(r)+∇⋅2Dn^t(r)ξ^t(r)(1)gdzie n^t(r) jest fluktuującą lokalną gęstością, która jest wielkością niekwilibrową: Gęstość lokalna cDFT jest jej uśrednioną fluktuacyjnie wartością dla układu równowagowego. Współczynnik D jest współczynnikiem dyfuzji w granicy małej gęstości (tzn. gdy w kąpieli znajduje się tylko jedna duża cząsteczka podlegająca ruchowi Browna), który może być obliczony na podstawie właściwości cząsteczek tworzących układ. Funkcję swobodnej energii F przyjmuje się jako funkcję Helmholtza z cDFT. Na poziomie hydrodynamiki (przed przyjęciem granicy nadtłumionej) jest ona bardziej ogólnie związana z gradientem ciśnienia w równaniach FH, a jej użycie jest rodzajem lokalnego przybliżenia równowagi powszechnego w niekwilibracyjnej mechanice statystycznej (26). Wreszcie, ξ^t(r) jest lokalnym białym szumem (z korelacjami funkcji delta w przestrzeni i czasie) wynikającym z małych cząsteczek kąpieli zderzających się z dużymi cząsteczkami i jest źródłem fluktuacji w modelu. Ważn± wła¶ciwo¶ci± tego modelu jest to, że zachowuje on liczbę cz±steczek przez cały czas, z wyj±tkiem granic układu. Tutaj koncentruję się na tym prostym, ale realistycznym modelu i pozostawiam podobny rozwój przy mniej restrykcyjnych założeniach do przyszłej pracy.

Użycie funkcji swobodnej energii cDFT w modelach stochastycznych zostało zakwestionowane (30), ponieważ to co występuje w typowych pochodnych modeli stochastycznych jest gruboziarnistą swobodną energią, a nie równowagową funkcją cDFT. Różnica wynika z fluktuacji wyraźnie reprezentowanych w modelu stochastycznym przez fluktuującą siłę i oczekuje się, że funkcja cDFT wynika z fluktuacyjnej średniej z gruboziarnistej funkcji. Tutaj, jak w prawie wszystkich zastosowaniach, funkcja swobodnej energii jest sumą złożonej funkcji twardej sfery i średniopolowego traktowania przyciągającego ogona potencjału (22). Można się spodziewać, że dla wkładów twardosferowych taka średnia będzie miała niewielki wpływ, ponieważ wszystkie korelacje są krótkozasięgowe. W przypadku układów z długozasięgowymi atrakcyjnymi ogonami można oczekiwać, że opis pola średniego stosowany tutaj i we wszystkich podobnych zastosowaniach będzie bardziej uzasadniony dla modeli gruboziarnistych niż dla modelu uśrednionego fluktuacyjnie, ponieważ to właśnie uśrednianie jest fizycznym źródłem efektów renormalizacji, które unieważniają opis pola średniego, np. w zjawiskach krytycznych. Stąd, kombinacja cDFT funkcji energii swobodnej twardej sfery i ogona przyciągającego pola średniego może być argumentowana jako dobre przypuszczenie co do gruboziarnistej energii swobodnej i złe przypuszczenie co do funkcji cDFT, a nie odwrotnie.

Możliwe byłoby wybranie jednej z funkcji modelu cDFT do użycia w modelu stochastycznym i przeprowadzenie bezpośrednich symulacji numerycznych. To podejście zostało ostatnio wykorzystane i jest obiecującą gruboziarnistą metodą symulacyjną do badania zarodkowania (31). Inną możliwością jest dalsze uszczegółowienie poprzez wprowadzenie zmiennych kolektywnych lub parametrów porządkowych. Ta droga była badana w innych miejscach, gdzie pokazano, że można odzyskać klasyczną teorię nukleacji (CNT) z odpowiednimi przybliżeniami (26, 32). Tak więc, teoria omawiana tutaj nie jest alternatywą dla CNT, ale raczej bardziej fundamentalną teorią, dla której CNT jest przybliżeniem. Tutaj jednak celem jest kontynuacja rozwoju teoretycznego przy użyciu narzędzi z teorii procesów stochastycznych i skupienie się na ścieżce zarodkowania jako podstawowym obiekcie.

Jeśli układ zaczyna się jako słaby roztwór (tj. stan podobny do pary) i spontanicznie zarodkuje klaster (albo kropla gęstego roztworu albo klaster krystaliczny), to początkowa lokalna koncentracja jest stała w całym układzie. Gdy klaster jest obecny, stężenie jest wysokie wewnątrz kropli, a niskie na zewnątrz: Różnica pomiędzy tymi stanami może być scharakteryzowana całkowicie w kategoriach lokalnej gęstości. Dla modeli stochastycznych, takich jak ten tutaj, możliwe jest podanie dokładnego wyrażenia na prawdopodobieństwo podążania dowolną ścieżką od pewnego początkowego rozkładu gęstości, n0(r), do dowolnego końcowego, nT(r), a poprzez poszukiwanie ścieżki z maksymalnym prawdopodobieństwem, można określić najbardziej prawdopodobną ścieżkę (MLP), która charakteryzuje przejście. Podstawowa idea pochodzi od Onsagera i Machlupa (33), a uogólnienie na dowolne procesy dyfuzyjne zostało podane przez Grahama (34). Jeśli stanem początkowym jest jednorodna faza macierzysta, a stan końcowy zawiera krytyczne (lub pokrytyczne) skupisko nowej fazy, to MLP będzie najbardziej prawdopodobną ścieżką nukleacji. W ogólności, wyznaczenie MLP jest wysoce nietrywialne, ale istotne uproszczenia mają miejsce w granicy słabego szumu (odpowiadającej fizycznie np. niższym temperaturom), w którym to przypadku teoria stochastyczna jest równoważna teorii dużych odchyleń Wentzella-Freidlina (35). Następnie można udowodnić (26), że MLP dla nukleacji musi przechodzić przez klaster krytyczny – fakt, który nie jest ogólnie prawdziwy w granicy silnego szumu i pokazuje, że zwykły obraz nukleacji ma zastosowanie tylko w tej granicy. Co więcej, można pokazać, że MLP może być skonstruowana przez rozpoczęcie w klastrze krytycznym i lekkie perturbacje układu w kierunkach niestabilnych tak, że deterministyczna część dynamiki∂∂tn^t(r)=D∇⋅n^t(r)∇δFδn^t(r)(2)powoduje spadek gęstości w dół gradientu energii swobodnej do fazy początkowej lub końcowej, w zależności od kierunku perturbacji. Złożenie tych dwóch cząstkowych ścieżek razem daje kompletne MLP. Zauważmy jednak, że mają one dwa bardzo różne charaktery: W rzeczywistości układ zaczyna się w fazie początkowej, a następnie jest napędzany przez fluktuacje w górę bariery swobodnej energii, aż do osiągnięcia klastra krytycznego, po którym następnie kontynuuje wzrost, aż do włączenia jak największej ilości materiału do nowej fazy. Druga część, rozpoczynająca się w klastrze krytycznym i rosnąca, jest tylko normalnym termodynamicznym wzrostem napędzanym przez gradient swobodnej energii. Pierwsza część, jednakże, jest napędzana fluktuacjami wbrew gradientowi swobodnej energii i jest wysoce nietrywialnym i użytecznym wynikiem, że MLP dla tego procesu może być określona przez opadanie „wstecz” w dół gradientu (26).

Zamiast bezpośrednio używać opadania gradientowego, tj. Eq. 2, od klastrów krytycznych, niniejsza praca wykorzystuje metodę łańcuchową (36), która jest matematycznie równoważna, ale oferuje zalety wydajności i prostoty. W szczególności, określa ona całą ścieżkę naraz i jest szczególnie użyteczna, gdy istnieje wiele pośrednich minimów swobodnej energii lub gdy gradienty swobodnej energii są słabe – oba te czynniki okażą się istotne poniżej. Szczegóły implementacji omówione są w Tekście Uzupełniającym, a tutaj zauważymy tylko, że w metodzie łańcuchowej pracuje się ze zbiorem gęstości, lub „obrazów”, rozmieszczonych wzdłuż całej ścieżki, przybliżając ją w ten sposób jako zbiór dyskretnych punktów. Przesuwając je zgodnie z równaniem 2 pod warunkiem zachowania równych odstępów między punktami, wyznaczana jest ścieżka opadania gradientowego. Metoda ta wymaga wstępnego odgadnięcia ścieżki, do czego wykorzystano prostą interpolację liniową pomiędzy punktami końcowymi. W niniejszych obliczeniach punktem wyjścia, czyli obrazem początkowym, jest jednolity układ o małej gęstości, a punktem końcowym na ścieżce jest klaster krytyczny. Obliczenia ujawniają, w jaki sposób układ ewoluuje od pierwszego do drugiego.

Wypada wspomnieć, że poszczególne elementy tej teorii były już wcześniej omawiane w podobnych kontekstach. Na przykład, w pracy Lutsko (24), bardzo podobna metoda nudged elastic band została użyta wraz z nowoczesną funkcją swobodnej energii cDFT do opisu nukleacji ciecz-ciecz. Jednakże, w tamtej pracy nie zrozumiano znaczenia wprowadzenia realistycznego opisu dynamicznego. Podobnie, badania pola fazowego, takie jak Qiu i Qian (37) oraz Backofen i Voigt (38, 39) wykorzystują prostsze funkcje swobodnej energii i dynamikę ad hoc wraz z metodą strunową. Podczas gdy krytyczne klastry są poprawnie wyznaczone, fizyczność ścieżek jest niejasna z powodu abstrakcyjnej natury dynamiki (np. brak lokalnego zachowania masy, gdy parametr porządku ma być interpretowany jako gęstość). Ponadto, efekty objętości wyłączonej, które dominują w strukturze na poziomie molekularnym, są poza domeną tych modeli. Obecna praca wykorzystuje zaawansowany model teorii miary fundamentalnej dla wkładu twardej sfery do funkcji swobodnej energii, który jak wiadomo daje bardzo dokładny opis układów twardej sfery w szczególności, a w połączeniu z modelem pola średniego dla przyciągającej części potencjału, struktury na poziomie molekularnym dla bardziej ogólnych potencjałów (22).

Podejścia do krystalizacji bardzo zbliżone w duchu do przedstawionego tutaj są również badane od pewnego czasu w społeczności PFC (25). Mogą one być rozumiane jako modele cDFT, które są uproszczone przez rozszerzenie funkcji swobodnej energii F wokół jednolitego stanu, n0, co skutkuje dwoma rodzajami terminów (41, 42). Pierwszy zestaw członów ma postać rozwinięcia gradientowego, które jest obcięte przy czwartym rzędzie. Drugi zestaw terminów ma postać rozwinięcia w zmiennej φ(r) ≡ (n(r) – n0)/n0, chociaż czasami unika się tego rozwinięcia (43). Oba rozwinięcia są niekontrolowane w stanie stałym, gdzie zmiany gęstości rzędu wielkości występują na najmniejszych skalach długości molekularnej i z tego powodu modele PFC nie mogą być parametryzowane do konkretnego potencjału oddziaływania. Z tego powodu, modele PFC nie mogą być sparametryzowane do konkretnego potencjału oddziaływania. Ta uproszczona funkcja swobodnej energii została sprzężona z dynamiką w wielu badaniach, w tym w niektórych, w których kombinacja jest bardzo podobna do tej zastosowanej tutaj (45), jako podstawa do zgrubnego klasyfikowania propagacji frontów krystalizacji w dwóch wymiarach (43) i nukleacji (44-46). Celem niniejszej pracy jest przezwyciężenie ograniczeń tych uproszczonych modeli przy jednoczesnym skupieniu się na ścieżce nukleacji poprzez koncepcję MLP.

.