Luận văn Nghiên cứu tính chất từ của vật liệu có cấu trúc nano trong điều kiện phân bố nhiệt không đều

ừ h ó a , M 11 1.4.2.1. Hiệu ứng GMR trong các màng mỏng đa lớp. Cấu trúc từ đa lớp (màng mỏng đa lớp) được tạo ra bằng phương pháp phóng xạ chùm ion, để tạo ra những lớp nguyên tử từ và không từ xen kẽ nhau. Vào năm 1988, nhóm của ALbert Fert của trường Đại học Tổng hợp Nam Paris đã nghiên cứu siêu mạng từ (001)Fe/(001)Cr chế tạo theo phương pháp epitaxy bằng chùm phân tử (MBE), các số liệu quan sát được về tính dẫn điện đã gây ấn tượng rất mạnh. Điện trở suất của hệ bão hoà khi ở trong từ tr- ường khoảng 20kg, giảm đi cỡ khoảng 2 lần so với khi không có từ trường ngoài đặt vào. Đây là một tỷ số MR lớn chưa từng có trước đó vì vậy hiệu ứng này đã được gọi là từ-điện trở khổng lồ. Có thể nói rằng, GMR xuất hiện từ sự tán xạ bởi tạp chất của điện tử dẫn nhưng với một số tính chất quan trọng là: có sự bất đối xứng spin, nghĩa là xác suất tán xạ của các spin - và spin- là khác nhau; và có sự tương quan về phương giữa spin và từ độ: song song hay phản song song. Khi chưa có từ trường tác dụng, từ độ trong các lớp xếp phản song song với nhau, hệ ở trạng thái điện trở suất cao. Khi có từ trường ngoài tác dụng, các véc tơ từ độ xoay lại theo hướng song song với nhau [12]. 1.4.2.2. Vai trò kích thước bề mặt. Một vấn đề đã được khẳng định là vai trò quan trọng đối với GMR của cấu trúc bề mặt phân cách trong các màng đa lớp. Nhiều nghiên cứu cho rằng độ nhám, hay độ mất trật tự, của bề mặt phân cách có vai trò làm tăng cường sự tán xạ phụ thuộc spin ở các bề mặt và do đó làm tăng biên độ của GMR. Những nghiên cứu khác lại cho rằng GMR có liên quan đến cấu trúc hoàn hảo của các lớp, những bề mặt trơn tru mới làm tăng độ lớn của GMR. Vì thế những bề mặt gồ ghề sẽ tạo ra những cầu nối ngẫu nhiên giữa các lớp từ làm tăng sự liên kết theo kiểu sắt từ. Do đó GMR giảm đi. Một ảnh hưởng lớn về vai trò của tạp chất trong các màng đa lớp là việc pha tạp có chọn lọc ở trong lớp từ hoặc xen giữa các bề mặt từ - không từ có thể đảo ngược hành vi của GMR, đó là hiệu ứng GMR ngược. Cũng trong nghiên cứu của Parkin vào năm 1990 đã nói ở trên, khi khảo sát một khoảng rộng các chiều dày khác nhau của các lớp trong hệ đa lớp Co/Cu, đã phát hiện ra hành vi biến thiên của GMR có tính chất với biên độ giảm dần khi tăng 12 chiều dày lớp Cu. Tính chất này đã được khẳng định trong nhiều hệ màng đa lớp khác khi thay đổi chiều dày lớp cách không từ. Sự thay đổi này của GMR đã phản ánh một hiện tượng xẩy ra giữa các lớp từ tương tự với tương tác trao đổi gián tiếp RKKY (Ruderman-Kittel- Kasuya- Yosida) giữa các nguyên tử tạp chất trong một đám các nguyên tử kim loại phi từ [13] . 1.4.2.3. Phương pháp đo GMR. Thường điện trở của của màng mỏng vẫn được tính theo công thức quen thuộc: 𝑅 = 𝜌 𝑙 𝑆 trong đó  là điện trở suất của vật liệu màng mỏng, 𝑙 là khoảng cách giữa hai điện cực đo điện áp và S là thiết diện của màng mỏng. Tuy nhiên, cần phải nhớ rằng, điện trở suất của vật liệu ở dạng màng mỏng khác với điện trở suất ở dạng khối do hiệu ứng giảm kích thước. Điện trở suất của màng mỏng có thêm thành phần do tán xạ bề mặt, 𝜌𝑠, ngoài ra các thành phần theo qui luật Matthiesen, gồm thành phần tán xạ bởi dao động mạng ( tán xạ phonon ), 𝜌𝑃𝐻, và do tạp chất cùng các loại sai hỏng khác, i. Do đó điện trở suất của màng mỏng được biểu diễn như sau:  = PH +i +i. Khi đó, điện trở suất của màng mỏng sẽ phụ thuộc vào chiều dày t của màng. Dạng phụ thuộc cụ thể giữa điện trở suất và chiều dày màng phụ thuộc vào tỷ số t/, với  là quãng đường tự do trung bình của điện tử. Cách đo điện trở thông dụng đối với trường hợp của màng mỏng là điện trở mặt, hay còn gọi là điện trở vuông, mà đơn vị tính bằng “ohm/ô vuông” đó là điện trở tính cho một ô vuông bất kỳ (không quan tâm đến kích thước) trên bề mặt màng mỏng. R=RS=/t() (1.4) Với t là chiều dày của màng mỏng RS là ký hiệu của điện trở vuông, không phụ thuộc vào các kích thước khác ngoài chiều dày của màng mỏng. Trong các công thức trên RS chính là điện trở vuông được định nghĩa ở trên. 1.4.2.4. Ứng dụng của spin + Các tính chất spin của điện tử có thể được dùng trong các máy tính lượng tử và thông tin lượng tử trong tương lai. + Trong thực sự định hướng của spin điện tử được sử dụng trong các cảm biến từ, đặc biệt là trong các đầu đọc và ổ cứng từ. 13 + Trong cơ học lượng tử và vật lý hạt, spin là một dạng nội tại của động lượng góc mang theo các hạt cơ bản, hạt tổng hợp (hadron) và hạt nhân nguyên tử. Spin là một trong hai loại động lượng góc trong cơ học lượng tử, loại còn lại là động lượng góc quỹ đạo. Toán tử xung góc quỹ đạo là đối trọng cơ học lượng tử với động lượng góc cổ điển của cuộc cách mạng quỹ đạo và xuất hiện khi có cấu trúc định kỳ đối với hàm sóng của nó khi góc thay đổi. Sự tồn tại của động lượng góc quay được suy ra từ các thí nghiệm, chẳng hạn như thí nghiệm Sternifer Gerlach, trong đó các nguyên tử bạc được quan sát thấy có hai mô men góc rời rạc có thể mặc dù không có động lượng góc quỹ đạo [14,15]. 1.5. ĐẤT HIẾM. Các nguyên tố đất hiếm và các kim loại đất hiếm, là tập hợp của mười bảy nguyên tố hóa học thuộc bảng tuần hoàn của Mendeleev, trái ngược với tên gọi (loại trừ promethi), có hàm lượng lớn trong Trái Đất. Các nguyên tố đất hiếm ở trong các lớp trầm tích, các mỏ quặng và cát đen. Nhóm đất hiếm thường không có tên trong sự sắp xếp khoa học. Tuy vậy, đất hiếm vẫn được tổ chức USPTO sắp xếp vào dạng hợp kim và các hợp chất khác, chính xác là nam châm đất hiếm từ các dạng khác nhau của nam châm. Các nguyên tố đất hiếm là : Xeri (Ce), dysprosi (Dy), erbi (Er), europi (Eu), gadolini (Gd), holmi ( Ho), lantan (La), luteti (Lu), neodymi (Nd), praseodymi (Pr), promethi (Pm), samari (Sm), scandi (Sc), terbium (Tb), thuli (Tm), ytterbi (Yb) và yttri (Y). Đất hiếm được gọi là "vitamin của nền công nghiệp hiện đại" nó là thành phần chủ yếu trong sản xuất các loại thiết bị và linh kiện trong công nghệ thông tin, y khoa, giao thông, hóa lọc dầu, luyện kim, quân sự và nhiều lĩnh vực khác. Đất hiếm có tính chất từ rất mạnh nên người ta thường sử dụng đất hiếm làm nam châm điện. Khả năng dẫn điện của đất hiếm rất tốt nên nó được sử dụng làm vật liệu siêu dẫn. Ngoài ra một số loại đất hiếm còn có tính phóng xạ mạnh nên nguy hiểm cho người sử dụng. Tóm lại, đất hiếm rất quan trọng vì không có chúng thì không thể sản xuất những sản phẩm thiết yếu phục vụ đời sống con người, và không thể thiếu trong chế tạo vũ khí, khí tài quân sự. Ứng dụng của đất hiếm. 14  Dùng để chế tạo các nam châm vĩnh cửu cho các máy phát điện  Dùng để đưa vào các chế phẩm phân bón vi lượng nhằm tăng năng suất và chống chịu sâu bệnh cho cây trồng  Dùng để chế tạo các nam châm trong các máy tuyển từ trong công nghệ tuyển khoáng  Dùng để diệt mối mọt, các cây mục nhằm bảo tồn các di tích lịch sử  Dùng chế tạo các đèn cathode trong các máy vô tuyến truyền hình  Dùng làm xúc tác trong công nghệ lọc hóa dầu và xử lý môi trường  Dùng làm vật liệu siêu dẫn  Các ion đất hiếm cũng được sử dụng như các vật liệu phát quang trong các ứng dụng quang điện  Được ứng dụng trong công nghệ laser [16]. 15 CHƯƠNG 2. CÁC MÔ HÌNH SPIN VÀ PHƯƠNG PHÁP MÔ PHỎNG MONTE-CARLO. 2.1 . CÁC MÔ HÌNH SPIN. Spin là một đặc trưng nội tại của hạt, nó giống như khối lượng và điện tích và được dùng để đặc trưng cho hạt đó. Một mô hình spin tổng quát bao gồm: i) Cấu trúc mạng tinh thể Cấu trúc mạng tinh thể biểu diễn trật tự sắp xếp hình học của các spin như mạng tinh thể 1 chiều, mạng tinh thể 2 chiều (mạng tam giác, mạng tổ ong ), mạng tinh thể 3 chiều (mạng lập phương, mạng lập phương tâm mặt). Căn cứ vào cấu trúc mạng, ta xác định được số spin tương tác với lân cận gần nhất (nearest neighbour), hay tương tác với next-nearest neighbor tại các nút mạng của tinh thể trong một số ví dụ dưới đây. Mạng spin một chiều mỗi spin sẽ tương tác với 2 spin lân cận gần nhất như trong hình vẽ 2.1 Hình 2.1 Mạng tinh thể 1D Mạng spin vuông 2 chiều, một spin sẽ tương tác với 4 spin lân cận gần nhất như trong hình vẽ 2.2. Hình 2.2 Mạng tinh thể vuông 2D Mạng tổ ong 2 chiều, mỗi spin sẽ tương tác với 3 spin lân cận gần nhất như trong hình vẽ 2.3 16 Hình 2.3 Mạng tổ ong 2 chiều. ii) Tập hợp các các giá trị của spin. Thông thường, spin có độ lớn 1S  , tùy thuộc vào trạng thái của spin trong mô hình, 1S  cho spin  và 1S   cho spin  trong mô hình spin Ising. Trong mô hình spin q -state Potts, spin nhận q giá trị khác nhau ( 2q  ). iii) Một định nghĩa về năng lượng tương tác giữa các spin. Trong vật liệu từ, tương tác trao đổi duy trì trật tự từ. Với mỗi mô hình spin, chúng ta cần đưa ra các giá trị các hằng số tương tác lân cận gần nhất hay hằng số tương tác kế tiếp lân cận gần nhất của các spin trong mạng tinh thể. iv) Một định nghĩa về tổng năng lượng của mô hình spin. Trình bày Hamiltonian bao gồm tổng năng lượng của hệ bao gồm năng lượng do các tương tác trao đổi giữa các spin và tương tác khi đặt hệ trong trường ngoài (từ trường). 2.1.1. Mô hình spin n -vector Mô hình spin n -vector hay mô hình ( )O n là một mô hình hết sức đơn giản trong số rất nhiều mô hình thuộc các chuyên ngành Vật lý, đặc biệt là trong Vật lý thống kê. Trong mô hình này, mỗi vector có độ dài bằng đơn vị, các spin iS đặt trên các nút mạng, Hamiltonian cho mô hình này là:   , .i j i j H J S S   (2.1) Trong đó tổng lấy theo tất cả các cặp spin lân cận ,i j và .i jS S biểu thị tích vô hướng Euclide thông thường. Các trường hợp đặc biệt của mô hình này là mô hình Ising ( 1n  ); mô hình XY ( 2n  ), mô hình Heisenberg ( 3n  ) 17 2.1.2. Mô hình spin Ising Mô hình Ising là một mô hình toán học về sắt từ trong cơ học thống kê. Mô hình gồm các spin có thể ở một trong hai trạng thái  và  . Xét một mạng tinh thể, ứng với mỗi vị trí i của mạng có một spin iS nhận giá trị 1 hoặc 1 . Hamiltonian cho mô hình spin này là , . .i j i j j j i j j H J S S h S      (2.2) Trong đó tổng đầu tiên chỉ lấy một lần cho mỗi cặp spin, jh là trường ngoài [16]. 2.1.3. Mô hình spin XY Mô hình spin XY cổ điển là một trong nhiều mô hình đơn giản trong cơ học thống kê, là trường hợp 2n  của mô hình n-vetor. Các spin cổ điển hai chiều iS nằm trên cùng một mặt phẳng. Spin tại các vị trí mạng x được biểu diễn như một vector trên một mặt hai chiều. Spin này được tham số hóa trong tọa độ cực, với giả thiết spin có giá trị đơn vị thì :  cos ( ),sin ( )S x x  (2.3) Với ( )x là định hướng của spin tại vị trí mạng x . Spin có thể hướng theo bất kỳ hướng nào, tương ứng với một số giá trị của  giữa 0 và 2 . Hamiltonian của mô hình spin XY với mô tả nói trên là :   , , . cosi j i j i j i j H J S S J        (2.4) Trong đó pha spin thứ i của i được đo từ trục ngang theo chiều kim đồng hồ và tổng lấy theo các spin lân cận. 18 2.1.4. Mô hình spin Heisenberg Các mô hình Heisenberg là trường hợp 3n  của mô hình n-vectơ. Có thể mô tả như sau: Lấy một mạng n chiều, và một tập các spin có độ dài đơn vị và mỗi spin được đặt tại một nút mạng. Mô hình được định nghĩa bới Hamiltonian:  , , .i j i j i j H J S S  (2.5) Với ,i jJ J nếu ,i j là các lân cận gần nhất, 0J  với mọi tương tác khác [17]. 2.2. PHƯƠNG PHÁP MÔ PHỎNG MONTE CARLO Mô phỏng Monte Carlo là phương pháp ngẫu nhiên để lấy mẫu trong một tập hợp thống kê. Phương pháp mô phỏng Monte Carlo cho phép mô hình hóa các hệ nhiều hạt có cấu trúc mạng tinh thể nhằm tính toán giá trị trung bình thống kê của các đại lượng đặc trưng cho hệ. Vật lý thống kê thường giải quyết các bài toán của hệ với nhiều bậc tự do. Vấn đề cơ bản trong vật lý thống kê là tính giá trị của những đại lượng vĩ mô quan sát được của hệ với giả thuyết đã biết là hàm Hamiltonia của hệ. Ví dụ, nếu một chất sắt từ có sự bất đẳng hướng theo một trục nào đó rất mạnh thì chúng ta có thể mô tả nó theo mô hình spin Ising, trong đó N spin iS tương tác theo Hamiltonian: , 1 i j i i j i H J S S h S      (2.6) Trong đó các spin iS tại nút mạng i có thể hướng lên  hoặc hướng xuống  dọc theo trục z , năng lượng trao đổi J trong công thức (2.1) chỉ được tính đối với các lân cận gần nhất và h là từ trường ngoài. Tuy nhiên, nếu chất sắt từ có sự bất đẳng hướng theo mặt phẳng (mô hình spin XY) hoặc bất đẳng hướng hoàn toàn (mô hình spin Heisenberg) thì Hamiltonian của hệ lúc này có dạng: 19   , x x y y x XY i j i j x i i j i H J S S S S h S     (2.7)   , . zHeisenberg i j z i i j i H J S S h S    (2.8) Từ mô hình cho bởi Hamiltonian, ta có thể tính toán được các giá trị của các đại lượng nhiệt động của hệ, ví dụ năng lượng trung bình E , hoặc độ từ hóa trung bình M được tính như sau. T H E N  (2.9) ii T S M N   (2.10) [18] 2.2.1. Thuật toán Metropolis Trong cổ điển, cấu hình của thuật toán Metropolis được tạo ra từ một trạng thái trước đó sử dụng xác suất chuyển dời cái mà phụ thuộc vào độ lệch năng lượng giữa trạng thái ban đầu và trạng thái cuối. Một chuỗi các trạng thái được tạo ra theo trật tự thời gian, nhưng thời gian trong mô phỏng Monte Carlo là thời gian Monte Carlo và không biết trước. Thuật toán Metropolis được sử dụng khá phổ biến trong mô phỏng Monte Carlo để tính giá trị trung bình của các đại lượng tuân theo một phân bố thống kê nào đó. n n n A P A (2.11) Với mô hình spin Ising ngẫu nhiên, sự phụ thuộc theo thời gian được mô tả bởi phương trình   ( ) ( )W ( )Wn n n m m m n m n P t P t P t t         (2.12) 20 Trong đó ( )nP t là xác suất tìm thấy của hệ ở trạng thái n taị thời điểm t , và Wn m là tốc độ chuyển từ trạng thái n sang trạng thái m . Khi hệ đạt trạng thái cân bằng ( ) 0nP t t   và hai số hạng ở vế phải là bằng nhau. Khi đó ta có ( )W ( )Wn n m m m nP t P t  (2.13) Xác suất của trạng thái n xảy ra trong hệ cổ điển cho bởi /( ) /n BE k TnP t e Z  (2.14) Trong đó Z là tổng thống kê của hệ. Xác suất này thường không được biết chính xác bởi mẫu số, chúng ta có thể tránh điều này bằng việc tạo ra một chuỗi Markov các trạng thái, cụ thể là tạo ra một trạng thái mới bằng chính trạng thái trước nó. Nếu chúng ta tạo ra trạng thái thứ n từ trạng thái thứ m, xác suất tương quan là tỉ lệ của xác suất cá nhân và sự xóa bỏ mẫu số. Dẫn đến, chỉ có độ lệch năng lượng của hai trạng thái là cần thiết n mE E E   (2.15) Tất cả những tốc độ chuyển nào thỏa mãn sự cân bằng chi tiết đều được chấp nhận. Sự lựa chọn đầu tiên của tốc đồ dời thường được sử dụng trong vật lý thống kê là dạng của Metropolis Khi 0E  : /1 0W BE k T m n e    (2.16) Khi 0E  : 1 0Wm n     (2.17) Trong đó 0 là thời gian yêu cầu để đạt được một spin lật. Thuật toán Metropolis được mô tả theo các bước như sau: Bước 1: Chọn một trạng thái ban đầu ngẫu nhiên Bước 2: Lật spin Bước 3: Tính giá trị độ lệch năng lượng E 21 Bước 4: Tạo số ngẫu nhiên r thỏa mãn 0 1r  Bước 5: Nếu / BE k Tr e , spin lật và cấu hình được cập nhật Sau đó lại lặp lại các bước từ bước 3 cho một vị trí nút mạng mới. 2.2.2. Phương pháp biểu đồ (Histogram) Ý tưởng sử dụng phương pháp biểu đồ để tăng lượng thông tin thu được từ phương pháp mô phỏng Monte Carlo. Trong khoảng 40 năm trở lại đây, phương pháp này đã và đang được ứng dụng thành công để nghiên cứu những hiện tượng tới hạn. 2.2.2.1. Phương pháp biểu đồ đơn. Phương pháp Monte Carlo tiến hành ở nhiệt độ 𝑇 = 𝑇0 tạo nên các cấu hình của hệ với tần suất tỷ lệ với trọng số Boltzmann exp (−𝛽0𝐻) với 𝛽0 = 1/𝑘𝐵𝑇0, và H là Haminton của hệ đang nghiên cứu. Để đơn giản ta xét mô hình Ising 3 chiều lân cận nhất với Hamiltonnian : 𝐻𝐼𝑠𝑖𝑛𝑔 = −𝐽 ∑ 𝑆𝑖𝑆𝑗〈𝑖,𝑗〉 . (2.18) Với nhiệt độ mô phỏng 𝑇0 tương ứng với hằng số liên kết 𝐾0 = 𝐽 𝑘𝐵𝑇0 , trọng số Boltzmann có thể viết dưới dạng exp (𝐾0𝐸), ở đây E là năng lượng của hệ. Xác suất tìm thấy hệ với năng lượng E và độ từ hóa 𝑀 = ∑ 𝑆𝑖𝑖 không thứ nguyên có dạng 𝑃𝐾0(𝐸, 𝑀) = 1 𝑍(𝐾0) 𝑊(𝐸, 𝑀)exp [𝐾0𝐸]. (2.19) Trong đó 𝑊(𝐸, 𝑀) là số các cấu hình (mật độ trạng thái) với năng lượng E và độ từ hóa M, và 𝑍(𝐾0) là hàm chia của hệ. Bởi vì mô phảng tạo ra các cấu hình tương ứng với phân bố xác suất cân bằng, biểu đồ H(E,M) của E và M giữ trong suốt quá trình mô phỏng gần đúng với hàm phân bố xác suất cân bằng, càng trở nên chính xác nếu số lần mô phỏng tiến đến vô hạn. Với một mô phỏng có thời gian hữu hạn,biểu đồ sẽ có sai số thống kê, nhưng đại lượng H(E,M)/N, với N là số lần đo đã thực hiện, vẫn chính xác đối với xác suất 𝑃𝐾0(E,M) trên khoảng giá trị của E và M được tạo ra trong suốt quá trình mô phỏng. Từ đó ta có thể viết lại phương trình như sau: 𝐻(𝐸, 𝑀) = 𝑁 𝑍(𝐾0) �̃�(𝐸, 𝑀)exp (𝐾0𝐸). (2.20) 22 Với �̃�(𝐸, 𝑀) là một xấp xỉ của mật độ trạng thái thực W(E,M). Từ sự giống nhau về dạng của phương trình (2.19) và (2.20), dễ dàng nhận thấy rằng nếu như ta biết được hàm phân bố tại một giá trị nào đó của K là có thể xác định các đại lượng tại các giá trị khác của K. Để thấy rõ điều này, chúng ta viết lại hàm phân bố cho các giá trị bất kỳ của K có dạng giống như (2.19) 𝑃𝐾(𝐸, 𝑀) = 𝑁 𝑍(𝐾0) 𝑊(𝐸, 𝑀)exp (𝐾0𝐸). (2.21) Tiếp theo, để ý là vì chúng ta đã biết được hàm phân bố tại nhiệt độ 𝐾0 , từ biểu đồ H(E,M), chúng ta có thể viết ngược lại phương trình (2.19) để xác định �̃�(𝐸, 𝑀) �̃�(𝐸, 𝑀) = 𝑍(𝐾0) 𝑁 𝐻(𝐸, 𝑀)exp (−𝐾0𝐸). (2.22) Nếu bây giờ ta thay 𝑊(𝐸, 𝑀) trong (2.21) bằng biểu thức cho �̃�(𝐸, 𝑀) trog (2.22), rồi lấy chuẩn hóa hàm phân bố, ta tìm được biểu thức liên hệ giữa biểu đồ đo được tại 𝐾 = 𝐾0 và hàm phân bố tại K bất kỳ. 𝑃𝐾(𝐸, 𝑀) = 𝐻(𝐸,𝑀)exp (𝐸∆𝐾) ∑ 𝐻(𝐸,𝑀)exp (−𝐾0𝐸)𝐸,𝑀 . (2.23) Với ∆𝐾 = 𝐾 − 𝐾0 . Từ 𝑃𝐾(𝐸, 𝑀), giá trị trung bình của bất kỳ hàm nào của E và M,f(E,M) có thể được tính như hàm liên tục của K 〈𝑓(𝐸, 𝑀)〉 = ∑ 𝑓(𝐸, 𝑀)𝐸,𝑀 𝑃𝐾(𝐸, 𝑀). (2.24) Phương trình (2.23) và (2.24) được xem là các phương trình biểu đồ đơn. Ưu điểm của phương pháp biểu đồ là có thể tính các đại lượng nhiệt động như một hàm liên tục của nhiệt độ K. Từ đó ta dễ dàng xác định được các cực trị của chúng. Sử dụng (2.24) ta cũng có thể tính được đạo hàm bậc nhất hay các bậc cao hơn của các đại lượng nhiệt động theo nhiệt độ. Ví dụ, đạo hàm của 〈|𝑚|〉 theo K có dạng 𝛿〈|𝑚|〉 𝛿𝐾 = (〈|𝑚|〉𝐸) − 〈|𝑚|〉〈𝐸〉). (2.25) Và đạo hàm bậc hai của nó là đơn giản 𝛿2〈|𝑚|〉𝐾 𝛿𝐾2 = (〈|𝑚|𝐸2〉𝐾 − 〈|𝑚|〉𝐾〈𝐸 2〉𝑘) − 2〈𝐸〉𝐾 𝛿〈|𝑚|〉𝐾 𝛿𝐾 . (2.26) 23 Các đạo hàm bậc cao hơn cũng có thể tính nếu cần. Khi một hàm, như nhiệt dung riêng, đạt đến giá trị cực đại thì đạo hàm của nó theo nhiệt độ K là bằng 0. Việc xác định các đỉnh này có thể chuyển sang việc tìm nghiệm của một phương trình. Ví dụ, bằng việc sử dụng phương pháp Newton dưới dạng các hàm tương quan giao nhau với năng lượng E như là hàm liên tục của K đã được mô tả ở trên. Quá trình này có thể tự động tìm đỉnh của các hàm một cách nhanh chóng với độ chính xác cao. 2.2.2.2 Phương pháp biểu đồ kép. Hầu hết phương pháp mô phỏng Monte Carlo đều sử dụng thuật toán Metropolis để tạo cấu hình Spin mới ở nhiệt độ cho trước. Sau đó người ta tính toán trung bình nhiệt động của các đại lượng quan sát được. Trong sự phát triển đáng ngạc nhiên gần đây, Ferrenberg và Swendson nhận ra là phương pháp này không sử dụng đầy đủ những dữ liệu có sẵn trong hàm phân bố thực của trạng thái đã thử. Họ đã phát triển những biểu đồ đơn của phương pháp Monte Carlo để tận dụng những lợi thế này nhằm thu nhận thêm thông tin. Khó khăn trong phương pháp biểu đồ đơn thực ra là do phân bố xác suất năng lượng khá hẹp tại mọi nhiệt độ . Điều này đã được cải thiện bằng phương pháp biểu đồ kép cũng được đưa ra bởi Ferrenberg và Swendson, trong đó dữ liệu thu được tại nhiều điểm nhiệt độ khác nhau kết hợp lại giúp cho hàm phân bố chính xác hơn các đại lượng quan sát được trên khoảng rộng của nhiệt độ. Để nghiên cứu kỹ về phương pháp biểu đồ kép, chúng ta xét Hamintonian cho trường hợp tổng quát có dạng : 𝐻(𝜎) = 𝐻0(𝜎) + 𝐾𝑆(𝜎) (2.27) Trong đó𝑆(𝜎) là toán tử năng lượng E hoặc là toán tử độ từ hóa M được xác định trên cấu hình spin{𝜎𝑖} của hệ, và K là nhiệt độ. Hàm chia của hệ có thể được viết dưới dạng : Z(K) = ∑ exp[H(σ)] = ∑ W(s)[KS]Sσ (2.28) ở đây W(s) là hàm mật độ trạng thái. 24 Xét R lần mô phỏng bằng phương pháp Monte- Carlo. Chúng ta thực hiện 𝑛𝑡ℎ mô phỏng tại nhiệt độ 𝐾𝑛 rồi lưu lại dưới dạng biểu đồ. {𝑁𝑛(𝑆)} là tổng số các giá trị {𝑛𝑛}. Bây giờ ta tính gần đúng của hàm chia: 𝑧𝑛(K) = ∑ 𝑁𝑛(𝑆)𝑒 (𝐾−𝐾𝑛)𝑆 𝑆 , (2.29) Hàm chia gần đúng này và hàm chia thực liên hệ với nhau qua biểu thức : 𝑧𝑛(𝐾)̅̅ ̅̅ ̅̅ ̅̅ = 𝑛𝑛𝑍(𝐾)/𝑍(𝐾𝑛) (2.30) Ta lại có biểu thức tính năng lượng tự do 𝑍(𝐾) − 𝑍(𝐾𝑛) = 𝑙𝑛𝑧𝑛(𝐾)̅̅ ̅̅ ̅̅ ̅̅ − 𝑙𝑛𝑛𝑛. (2.31) Hàm mật độ trạng thái liên hệ đến biểu đồ: 𝑊(𝑆) = 𝑁𝑛(𝑆)/𝑛𝑛̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅̅ ̅𝑒 𝑓𝑛−𝐾𝑛𝑆 (2.32) Trong đó, 𝑓𝑛 = 𝐹(𝐾𝑛) là tham số bằng năng lượng tự do tại 𝐾𝑛, và nó sẽ được tính bằng phương pháp tự hợp. Chúng ta thực hiện mô phỏng trên tập các giá trị {𝐾𝑛|𝑛 = 1, 𝑅}, chúng ta có thể kết hợp chúng lại thành một biểu thức dưới dạng tổng quát cho W(S) W(S) = ∑ pn(S)Nn(S)nn −1efn−KnSRn=1 (2.33) Với ∑ pn(S) = 1 R n=1 (2.34) Nếu chúng ta thay thế các biểu đồ vào công thức (2.33) và cực tiểu sai số chúng ta tìm được : 𝑝𝑛𝑆 = 𝑛𝑛𝑒 𝐾𝑛𝑆−𝑓𝑛 ∑ 𝑛𝑚𝑒𝐾𝑚𝑆−𝑓𝑚 𝑅 𝑚=1 (2.35) Nếu ta định nghĩa 𝑃(𝑆, 𝐾) = 𝑊(𝑠)𝑒𝐾𝑆 Ta thu được phương trình biểu đồ kép 𝑃(𝑆, 𝐾) = ∑ 𝑁𝑛(𝑆)𝑒 𝐾𝑆𝑅 𝑛=1 ∑ 𝑛𝑚𝑒𝐾𝑚𝑆−𝑓𝑚 𝑅 𝑚=1 (2.36) ở đây 𝑒𝑓𝑛 = ∑ 𝑃(𝑆, 𝐾𝑛)𝑆 25 Giá trị trung bình của bất kỳ toán tử nào của S cũng có thể tính được như là một hàm liên tịc theo K: 〈𝐴(𝑆)〉(𝐾) = ∑ 𝐴(𝑆)𝑃(𝑆, 𝐾)/𝑧(𝐾)𝑆 (2.37) 26 Trong đó 𝑧(𝐾) = ∑ 𝑃(𝑆, 𝐾) 𝑆 Giá trị của 𝑓𝑛 được tính bằng phương pháp tự hợp của hai phương trình (2.36) và (2.37). Để tăng tốc độ hội tụ trong quá trình tính lặp, ta sử dụng đạo hàm của giá trị mới 𝑓𝑛 theo giá trị cũ của nó. 2.2.3. Phương pháp Wang - Landau. Như chúng ta đã biết hầu hết các thuật toán trong mô phỏng Monte Carlo đều được thực hiện tại một nhiệt độ xác định. Bởi vậy, các phương pháp trước đây, hoặc là kém chính xác (phương pháp cổ điển), hoặc là chỉ có hiệu quả cao trong việc tính toán các hiệu ứng kích thước hữu hạn (phương pháp biểu đồ đơn và kép). Các phương pháp này được ứng dụng rộng rãi và khá hiệu quả đối với chuyển pha loại II, nhưng lại rất khó khăn khi nghiên cứu chuyển pha loại I. Việc xác định đỉnh kép trong hàm phân bố năng lượng cũng như xác định độ rộng tatent giữa hai đỉnh kép đó đòi hỏi phải thực hiện mô phỏng tại nhiệt độ rất gần điểm chuyển pha, tức là phải tiến hành mô phỏng tại rất nhiều điểm nhiệt độ quanh điểm chuyển pha với độ chính xác cao. Mới đây, Wang và Landau đã đề xuất một thuật toán mới gọi là thuật toán biểu đồ “phẳng”, còn gọi là phương pháp Wang - Landau. Khác với các thuật toán mô phỏng Monte Carlo trước đây, phương pháp Wang -Landau mô phỏng trên không gian năng lượng nhằm xác định hàm mật độ trạng thái. Từ đó, ta có thể xác định các định các đại lượng vật lý như là hàm liên tục theo nhiệt độ tương tự như trong phương pháp giản đồ đơn (kép). Tiếp theo đây, chúng ta có thể tìm hiểu chi tiết thuật toán này. Thuật toán này dựa trên quan sát thấy rằng, nếu như ta thực hiện một bước đi ngẫu nhiên trong không gian năng lượng với xác suất tỷ lệ nghịch với nghịch đảo mật độ trạng thái 1/g(E), lúc đó biểu đồ phẳng được tạo ra từ phân bố năng lượng. Quá trình này được thực hiện bằng cách điều chỉnh mật độ trạng thái ước định theo theo phương pháp có hệ thống, nhằm thu được biểu đồ phẳng trên toàn bộ vùng năng lượng cho phép và đồng thời làm cho mật độ trạng thái hội tụ tiến dần đến giá trị đúng của nó. Ban đầu, giá trị của mật độ trạng thái chưa được xác định, vì vậy ta có thể đặt g(E)=1 cho tất cả các năng 27 lượng E. Sau khi đặt giá trị ban đầu cho mật độ trạng thái, ta thực hiện mô phỏng theo bước đi ngẫu nhiên. Trong trường hợp tổng quát , nếu như năng lượng trước và sau khi lật spin tương ứng là 𝐸1 sang 𝐸2 có dạng : 𝑃(𝐸1→𝐸2) = min ( 𝑔(𝐸1) 𝑔(𝐸2) ). (2.38) Đây cũng chính là xác suất lật spin. Mỗi lần ta tìm thấy trạng thái của hệ có năng lượng E, ta làm mới mật độ trạng thái ứng với năng lượng đó bằng cách nhân mật độ trạng thái bởi hệ số điều chỉnh f>1, nghĩa là : 𝑔(𝐸𝑖) → 𝑔(𝐸𝑖)𝑓. Hệ số điều chỉnh ban đầu được chọn khá lớn = 𝑓0 = 𝑒 1 = 2.71828 . Nó cho phép chúng ta tìm đến một cách nhanh c