Átlagtér számolás Mathematicaval

Radnóczi Katalin fizikus, II. Ph.D. hallgató

2004. január

Bevezetés

A szilárdtestfizikában gyakran használják az ún.másodkvantált formalizmust,amelynek lényege,hogy a lehetséges elektronállapotokat (amelyeket a kvantummechanika a hullámfüggvényekkel, mint egy Hilbert-tér elemeivel ír le) a vákuum állapotból (amely nem tartalmaz semmilyen részecskét sem) egy (a megfelelõ állapotba) keltõ operátorral hozzák létre. A különbözõ állapotokba keltõ operátorok és ezek adjungáltjai,az eltüntetõ operátorok egymással meghatározott antikommutációs relációkat elégítenek ki. Így biztosítják többek között, hogy két elektron ne kerülhessen azonos állapotba. A rendszerben minden fizikai mennyiség egy operátorral fejezhetõ ki,amely a fenti keltõ és eltüntetõ operátorokból áll elõ. Minket általában ezek sajátvektorrendszere érdekel. Például a rendszer alapállapota az energia (Hamilton) operátorának legkisebb sajátértékû sajátállapota.
Ebben a programban egy ilyen, nemkommutatív algebrán alapuló, Mathematicaban megírt programot fejlesztek tovább egy átlagtér számoláshoz.Ez az alábbi részekbõl épül fel:

Definiálom a fent említett nemkommutatív algebrát operátorokra

Olyan függvényt írok, amelyek egy operátorhoz hozzárendeli az adjungáltját, illetve egy állapot(vektor)hoz a duális tér megfelelõ elemét.

Definiálom az elektronállapotoknak azon leszûkített bázisát, amelyre szükségem van, azaz a Hilbert-teret

A már definiált operátorok szorzataként elõállítom egy speciális probléma Hamilton-operátorát. Ezzel kiszámolom a Hamilton operátor mátrixát egy konkrét bázison,ami a kölcsönhatási paraméterekkel megadva túl bonyolult lesz.

Olyan függvényt írok, ami egy operátorhoz az ún. átlagtér-szétcsatolt operátorát rendeli hozzá (ahol négy keltõ/eltüntetõ operátor szorzata szerepel bennük, ott kettõt a várható értékével helyettesítek. Ez azért egyszerûsíti a problémát, mert az -- új paraméterként bevezetett -- várható érték nem operátor, hanem csak egy szám).

Az elõzõ függvényt alkalmazom a Hamilton operátorra, és meghatározom a mátrixát az adott bázison

Ezt diagonalizálom, és meghatározom az állapotösszeget (a mátrix súlyozott spurját)

Felírom az önkonzisztencia-egyenletet a várható érték(ek)re (ebben az állapotösszeg deriváltja szerepel), melyet megoldok ezekre a paraméterekre,amelyek bizonyos fizikai mennyiségek várható értékeit közelítik

Nemkommutatív algebra

<< "/home/kati/fizika/nb/NonCommutativeMultiply.m";

ClearAll[delta];
delta[i1_Integer, i2_Integer] := If[i1 == i2, 1, 0];
delta[i1_List, i2_List] := 
        If[SameQ[i1, i2], 1, 0];

    
ClearAll[nce];
nce[x_] := Expand[NonCommutativeExpand[x]];

Két (elektronállapotok Hilbert terén ható) operátor kommutátora és antikommutátora

\(\(\( (*ClearAll[commute, \ 
        acommute]; 
commute[x_, y_] := 
        nce[x ** y - y ** x]; 

      acommute[x_, y_] := 
        nce[x ** y + 
            y ** x];*) \)\(
\)\(
\)\(\(ClearAll[comm];\)

    \(Options[
          comm] = {AntiCommutator -> True};\)
 (*comm[
          x_, y_, opts___?
            OptionQ] := \((kif = 
              nce[\((x)\) ** \((y)\) + \((y)\) ** \((x)\)]; 
            If[\(AntiCommutator /. {opts}\) /. Options[comm], kif, 
              MapAt[\(-1\) # &, 
                kif, {{2, 1}}]])\)\ itt\ az\ a\ probléma, \ 
      hogy\ az\ elõtt\ kellene\ a\ MapAt - 
        et\ végrehajtani, \ 
      mielõtt\ x - et\ és\ y - t\ behelyettesíti*) 

    \(comm[x_, y_, opts___?OptionQ] := 
        If[\(AntiCommutator /. {opts}\) /. Options[comm], 
          nce[x ** y + y ** x], nce[x ** y - y ** x]];\)\)\)\)

Elektronokat keltõ (fc[]) és eltüntetõ (fa[]) operátorok definiálása.  Annyi ilyen operátor van, ahány állapotban lehet egy elektron, ezért ezeket függvényként definiálom, melynek argumentuma a megfelelõ állapot. Általában végtelen sok ilyen állapot lehetséges, mi itt egy leszûkített, véges dimenziós Hilbert térrel, tehát véges fajta ilyen operátorral dolgozunk. 
Az ilyen állapotok valamelyikében levõ elektronokra ható bármilyen összetett operátort ezen alapvetõ operátorokból felépített szorzatokkal és ezek összegeivel meg lehet adni.
Teljesül, hogy egy azonos állapotba keltõ és eltüntetõ operátor egymás adjungáltjai.

ClearAll[fa, fc];
SetNonCommutative[fa, fc];

A kvantummechanika szerint két elektron nem lehet ugyanabban az állapotban (Pauli elv). Ezt másodkvantált formalizmusban úgy fejezzük ki, hogy egy olyan összetett operátor, amely két, azonos állapotba keltõ operátort tartalmaz, úgy hat egy tetszõleges elektronállapotra, hogy 0-ba viszi. Az ilyen operátorok tehát 0-val helyettesíthetõk. Ugyanez igaz két azonos állapotból eltüntetõ operátor szorzatára is, mert ez az elõzõ operátor adjungáltja.

fa /: z___ ** fa[x__] ** fa[x__] ** w___ := 0;
fc /: z___ ** fc[x__] ** fc[x__] ** w___ := 0;

A másodkvantált formalizmusban definiált keltõ és eltüntetõ operátorok nem felcserélhetõk egymással, hanem meghatározott antikommutációs relációkat elégítenek ki, amelyek az alábbiak:

fa /: 
        z___ ** fa[x__] ** fa[y__] ** 
          w___ := -z ** fa[y] ** fa[x] ** 
              w /; ! (OrderedQ[{{x}, {y}}]);
fc /: 
        z___ ** fc[x__] ** fc[y__] ** 
          w___ := -z ** fc[y] ** fc[x] ** 
              w /; ! (OrderedQ[{{x}, {y}}]);

     (*fc /: 
          z___ ** fa[x__] ** fc[y__] ** 
            w___ := -nce[
              z ** (fc[y] ** fa[x] - delta[{x}, {y}]) ** 
                w];*) 
fc /: comm[fa[x___], fc[y___]] := delta[{x}, {y}];

A Hilbert tér egy speciális eleme a vákuum állapot (vac), amely nem tartalmaz elektront, ebben nincsen tehát semmilyen részecske. A Hilbert tér duálisában létezik ennek a megfelelõje (vaca), amelyet a vákuum állapot adjungáljának is szoktak nevezni. 
Úgy definiáljuk, hogy 1-re normált legyen. Ez egyben azt is jelenti, hogy a belõle egy db. keltõ operátorral letrehozott egy elektront tartalmazó állapotok (pl. a "2." állapot, mint Hilbert-tér-elem fc[2]**vac alakban áll elõ) normája is 1. Több elektront tartalmazó állapotokat is létrehozhatunk, ha a vákuum állapotra több keltõ operátort is alkalmazunk (pl. fc[1]**fc[2]**vac egy olyan állapot, amelyben az 1. és a 2. állapotban is van egy-egy részecske).
A vákuum állapot kitüntetett tulajdonsága, hogy nem tartalmaz egyetlen részecskét sem, ezért nem lehet belõle eltüntetni elektronokat, tehát egy eltüntetõ operátorral hatva rá, 0-t kell, hogy kapjunk. (Ennek az "adjungáltja" is igaz: vaca-ra jobbról hatva egy keltõ operátorral 0-t kapunk.)
Ha egy állapotot balról megszorzunk egy másik állapot adjungáltjával, akkor a két állapot skalárszorzatát kapjuk, ami egy szám. A saját adjungáltjával való szorzat az állapot normájának négyzete. Ha jobbról szorozzuk meg egy másik állapot adjungáltjával, akkor egy (egyelektron-állapotok terén ható) operátort kapunk. 
Az fc[2]**vac állapot adjungáltja (a duális tér megfelelõ eleme): vaca**fa[2].

\(\(\( (*the\ vacuum\ \(state : \ 
          element\ of\ the\ Hilbert\ space\ containing\ no\ electrons\)*) \)\(
\)\(\(ClearAll[vac, 
        vaca];\)

    \(SetNonCommutative[vac, vaca];\)

    \(vac /: z___ ** fa[__] ** vac := 0;\)

    \(vaca /: vaca ** fc[__] ** z___ := 0;\)

    \(vaca /: vaca ** vac := 1;\)

    \(vaca /: z___ ** vaca ** vac ** w___ := z ** w;\)

    \(vaca /: z___ ** vaca ** fa[x___] ** fc[y___] ** vac ** w___ := 
        z ** delta[{x}, {y}] ** w;\)

    \(vaca /: vaca ** fa[x___] ** fc[y___] ** vac := 
        delta[{x}, {y}];\)\)\)\)

Az x állapotban levõ részecskék számának operátora (értéke a Pauli elv alapján 0 v. 1 lehet)  egyelektron-állapotokra:

ClearAll[n];
n[x__] := fc[x] ** vac ** vaca ** fa[x];

A delta fv. helyett nem hasznalható a KroneckerDelta, mert:

KroneckerDelta[{1, b}, {a, b}]

comm[fa[1], fc[1], AntiCommutator -> False]

comm[fa[1], fc[2]]

A definiált operátorok helyes mûködésének ellenõrzése

n[1]

nce[1]

vac ** vac

nce[fc[1] ** vac ** vaca ** fc[2]]

nce[fc[1] ** vac ** vaca ** fa[1]]

nce[vaca ** fa[1] ** fc[1] ** vac]

vaca ** fa[1] ** fc[1] ** vac

vaca ** vac

Operátorok adjungálása

Megfordítjuk az operátorok sorrendjét, és mindegyiket külön-külön adjungáljuk (keltõ-> eltüntetõ es fordítva). 
A komplex szám szorzófaktorokat konjugálni kell.

FullForm[(2 + 3  I) fa[1] ** fc[2]]

Apply[List, ((2 + 3  I) fa[1] ** fc[2])]

Ismeretlen összetett kifejezések adjungálása, amelyek komplex számokat tartalmaznak: a bennük szereplõ komplex számokat konjugáljuk

szamadj[y_] := y /. x_Complex -> Conjugate[x]

szamadj[(2 + 3  I) s + (4 + 5  I) 
          k]

Az opadj függvény csak olyan operátorokat tud adjungálni, amelyekben nem szerepel komplex szorzófaktor, és csak egy nemkommutatív operátorokból álló szorzat sorozatok (melyek egyelemûek is lehetnek). 
Az adj függvény már ilyenek tetszõleges komplex számokkal vett lineáris kombinációját is helyesen adjungálja.

ClearAll[opadj, adj];
opadj[x_] := 
        Module[{}, 
          If[Length[x] > 1, 
            NonCommutativeMultiply @@ (Reverse[
                    List @@ x] /. {fa :>  fc, 
                    fc :>  fa, vac :>  vaca, 
                    vaca :>  vac}), 
            x /. {fa :>  fc, fc :>  fa, 
                vac :>  vaca, 
                vaca :>  vac}]];
adj[z_] := 
      z /. {y_NonCommutativeMultiply :> opadj[y], 
          x_Complex -> Conjugate[x], fa :>  fc, 
          fc :>  fa, vac :>  vaca, 
          vaca :>  vac}

kif = 
      Complex[2, 3]*fc[3] ** fc[4] + 3  fc[1] ** vac + 
        2  I fc[2] ** fa[5] + vaca + 5  fc[7]

opadj[kif]

opadj[fa[1] + fc[2]]

szamadj[kif]

Csak az adj[] ad helyes eredményt tetszõleges összetett operátor esetén, mint amilyen kif

adj[kif]

Conjugate[kif]

opadj[fa[1] ** fc[2] ** vac]

adj[fa[1] ** fa[2] ** fa[3]]

adj[fa[1] ** fc[2]]

nce[fc[1] ** fa[2] ** adj[fc[1] ** fa[2]]]

adj[fa[1]]

A részecskeszám-operátor önadjungált:

adj[n[1]]

adj[n[3]] - n[3]

Az igy definiált adjungáló operátor állapotok adjungálására is alkalmas:

adj[fc[2] ** vac]

adj[adj[vaca ** fa[2]]]

adj[vac]

Skalárszorzatok és normanégyzetek számítása

adj[vac] ** vac

adj[fc[2] ** vac] ** fc[2] ** vac

adj[fc[2] ** vac] ** fc[3] ** vac

Bázisállapotok

Jelen esetben egy olyan problémát tekintek, ahol az elektronok egy egydimenziós lánc rácshelyein foglalhatnak helyet, minden rácshelyen nyolcféle állapotban lehetnek (α. vagy β. atomon, a vagy b pályán, fel(u) vagy le(d) spinállapotokban). A rácshely a keltõ/elt. operátorok elsõ argumentuma, a többi az állapotra vonatkozik azon a rácshelyen. Az egyelektron-állapotokat a vákuumállapotból egy darab elektronkeltõ operátorral hozzuk létre, az alábbi állapotban például egy elektron van az 1. rácshelyen, az α. atomon, a pályán, felfelé mutató (u) spinnel:

fc[1, α, a, u] ** vac

Egy elektron lehetséges állapotai a j. (bármely) rácshelyen (a j index jelen esetben egy egydimenziós lánc rácshelyein fut végig: j ∈ {1,2, ..., N}, ahol N a rácshelyek száma)

atomlist = {α, β};
orblist = {a, b};
spinlist = {u, d};

ClearAll[st];
st[j_] := 
      Partition[Flatten[Outer[List, {j}, atomlist, orblist, spinlist]], 
        4]

st[2]

Egy elektron lehetséges állapotai az 1. és a 2. rácshelyen

stfor2sites = 
      Partition[Flatten[Outer[List, {1, 2}, atomlist, orblist, spinlist]], 
        4]

Késõbbi matrixelem-számítashoz egy φ vektort hozunk létre, ami az 1. rácshely lehetséges állapotaiból áll, és egy másikat (φa), ami ezek adjungáltjait tartalmazza. O operátorok mátrixait ezen az altéren φa**O **φ alakban számolhatjuk majd ki. A vektorbeli minden állapot egyre normált.

φ = # ** vac & /@ Apply[fc, st[1], {1}]

φa = adj /@ φ

Egy elektron 1. és a 2. rácshelyen felvett lehetséges állapotainak vektora:

Inner[NonCommutativeMultiply,φa,φ,Plus]

Két elektron (az egyik az 1., a másik a 2. rácshelyen) felvett lehetséges állapotainak vektora (8*8=64 féle lehet):

st21 = 
        Partition[Flatten[Outer[List, {1}, atomlist, orblist, spinlist]], 
          4];
st22 = 
        Partition[Flatten[Outer[List, {2}, atomlist, orblist, spinlist]], 
          4];

st21

φ22 = -Flatten[
          Outer[NonCommutativeMultiply, 
            Apply[fc, st[2], {1}], φ]]

φ22a = adj /@ φ22;

Hamilton operátor

A rendszer energiaoperátora H=(∑_i)((Subscript[H, i, i + 1]) + (Subscript[E, i]))alakú, ahol i a rácshelyeken fut végig (i ∈ {1,2, ..., N}), (Subscript[H, i, i + 1])két szomszédos rácshelyet köt össze (a szomszédos rácshelyeken levõ elektronok kölcsönhatását fejezi ki, amely attól függ, hogy az egyes rácshelyeken levõ részecskék milyen állapotban vannak), (Subscript[E, i]) egy-rácshely operátor. Mindkét fajta operátor általában bonyolult alakú, különbözõ spin és/vagy pálya kölcsönhatásokból épül fel, melyek alakját a rács/spin/idõtükrözési szimmetriák szabályozzák. Hogy (Subscript[E, i])-t felírhassuk, elõször egyszerûbb, elemi operátorokat definiálunk.

Elemi operátorok

atomlist = {α, β};
orblist = {a, b};
spinlist = {u, d};

Thread[
      Apply[fc, 
          Partition[Flatten[Outer[List, {i}, atomlist, {a}, spinlist]], 
            4], {1}] ** vac ** vaca ** 
        Apply[fa, 
          Partition[Flatten[Outer[List, {i}, atomlist, {b}, spinlist]], 
            4], {1}]]

Elemi operátorok (atomokra/pályákra/spinekre)

ClearAll[a, b, leptetooper, xoper, yoper, 
      zoper]
leptetooper[{x___}, {y___}] := 
        Plus @@ Thread[
            Apply[fc, Partition[Flatten[Outer[List, {j}, x]], 4], {1}] ** 
              vac ** vaca ** 
              Apply[fa, 
                Partition[Flatten[Outer[List, {j}, y]], 
                  4], {1}]];
xoper[{x___}, {y___}] := (leptetooper[{x}, {y}] + 
              leptetooper[{y}, {x}])/2;
yoper[{x___}, {y___}] := (leptetooper[{x}, {y}] - 
              leptetooper[{y}, {x}])/(2  I);
zoper[{x___}, {y___}] := (leptetooper[{x}, {x}] - 
              leptetooper[{y}, {y}])/2;

A fent definiált léptetõ operátorok hatása az elektronállapotokra: kétféle állapot között "léptetnek", azaz az alábbi, pályákra vonatkozó operátor például egy "b" pályán levõ elektront egy "a" pályaállapotba visz át, miközben a spinjét, és azt, hogy az "α" vagy a "β" jelzésu atomon tartózkodik, változatlanul hagyja. A pályákra vonatkozó operátorokból három van, amelyek egymástól függetlenek, ezeket (Subscript[τ, x]),  (Subscript[τ, y])és  (Subscript[τ, z])jelöli. Ezek közül az elsõ kettõ az a->b, b->a operátorok bizonyos lin. kombinációi, míg a harmadik nem léptet, csak ±(1/2)-del szorozza az állapotot a szerint, hogy "a" vagy "b" pálya szerepel benne.
Hasonló operátorokat konstruálhatunk, melyek az atom, ill. spinállapotokat fogják manipulálni. A spinállapotok között (lefelé) léptetõ operátorok pl. az "u" állapotokat "d" állapotba viszik, és az atom/pálya állapotokat hagyják változatlanul. Ezekbõl  (Subscript[s, x]),  (Subscript[s, y])és  (Subscript[s, z])operátorokat építek fel, az elõzõekhez teljesen hasonlóan. 
Az atom-állapotokra vonatkozó, az elõzõek analógiájára felépített operátorok, (Subscript[η, x]),  (Subscript[η, y])és  (Subscript[η, z])pedig pl. α->β manipulációkat is tartalmaznak.

leptetooper[{atomlist, {a}, spinlist}, {atomlist, {b}, 
        spinlist}]

xoper[{atomlist, {a}, spinlist}, {atomlist, {b}, spinlist}]

Elemi operátorok (pályákra)

τx[i_] := 
        xoper[{atomlist, {a}, spinlist}, {atomlist, {b}, spinlist}] /. 
          j -> i;
τy[i_] := 
        yoper[{atomlist, {a}, spinlist}, {atomlist, {b}, spinlist}] /. 
          j -> i;
τz[i_] := 
        zoper[{atomlist, {a}, spinlist}, {atomlist, {b}, spinlist}] /. 
          j -> i;

τy[2]

Az így definiált 9 operátor önadjungált, pl.

Simplify[% - adj[%]]

Elemi operátorok (atomokra)

ηx[i_] := 
        xoper[{{α}, orblist, spinlist}, {{β}, orblist, 
              spinlist}] /. j -> i;
ηy[i_] := 
        yoper[{{α}, orblist, spinlist}, {{β}, orblist, 
              spinlist}] /. j -> i;
ηz[i_] := 
        zoper[{{α}, orblist, spinlist}, {{β}, orblist, 
              spinlist}] /. j -> i;

Elemi operátorok (spinekre)

sx[i_] := 
        xoper[{atomlist, orblist, {u}}, {atomlist, orblist, {d}}] /. 
          j -> i;
sy[i_] := 
        yoper[{atomlist, orblist, {u}}, {atomlist, orblist, {d}}] /. 
          j -> i;
sz[i_] := 
        zoper[{atomlist, orblist, {u}}, {atomlist, orblist, {d}}] /. 
          j -> i;

sz[1]

Ezen 9 elemi operátornak fizikai jelentése van, pl. a az i. rácshelyen levõ elektron spinje (az elektron saját impulzusmomentuma, vagyis perdülete) z irányú komponensének operátora. A kvantummechanika szerint úgy számolhatjuk ki, hogy bizonyos állapotokban mennyi a spin z komponensének várható értéke (melyet h/2π (h a Planck állandó) egységekben kell érteni), ha az operátort az állapottal jobbról, az állapot adjungáltjával balról megszorozzuk. például:

nce[
      adj[fc[1, α, b, u] ** vac] ** sz[1] ** fc[1, α, b, u] ** 
        vac]

A spinoperátor a kvantummechanikából ismert, és gyakran használatos. A τ és η operátorok jelentése nem ilyen egyértelmû és szembetünõ, a pályadegenerált rendszerek tárgyalásánál mégis gyakran bevezetik õket a spinoperátorok analógiájára. Az elektron spinje ugyanis kétféle lehet, és jelen esetben a pályaállapotok ill. az atomok száma, melyen az elektron tartózkodhat egy rácshelyen, szintén kettõ.

Ez a 9 elemi operátor egymással kommutál, például:

comm[τx[1], sx[1], AntiCommutator -> False] // 
      Simplify

Másodrendû szorzataikbol újabb 54 operátor epíthetõ fel (egy adott rácshelyen). A Hamilton operátor 2 rácshelyre ((Subscript[H, ij]))a szomszéd rácshelyeken levõ 64*64 fajta operátorszorzatókbol épül fel, amelyek közül azonban a szimmetria nem mindenfajta szorzatot enged meg.    
Egy rácshelyhez tartozó operátor ((Subscript[H, i])) például:

(τx[1]) ** (sx[1])

nce[(τx[1]) ** (sx[1])] // Simplify

nce[ηy[1] ** τx[1] ** sz[1]] // Simplify

Hamilton operátor egy rácshelyen

A Hamilton operátor az imént definiált elemi operátorokból és azok szorzataiból épül fel. 
A szimmetria szerinti osztályozás alapján a Hamilton operátorban elsõrendû tagok:

h[i_] := 
      nce[2  λ*τy[i] ** sz[i] + a2*ηx[i] + 
          a3*τy[i] ** sz[i] ** ηx[i] + 
          a4*(τz[i] ** sx[i] - τx[i] ** sy[i]) ** ηy[
                i]]

h[1] // Simplify

Ennek a sajátértékeit keressük az egy-rácshely bázison:

φ

hmat = 
        Table[nce[φa[[i]] ** h[1] ** φ[[j]]], {i, 
            8}, {j, 8}];

% // MatrixForm

Eigenvalues[hmat] // Simplify

Csak a furcsa tag sajátállapotai

Table[
      nce[φa[[i]] ** ((τz[1] ** sx[1] - τx[1] ** 
                    sy[1]) ** ηy[1]) ** φ[[j]]], {i, 
        8}, {j, 8}]

Eigensystem[%] // Simplify

Két rácshely közötti kölcsönhatást leíró Hamilton operátor átlagtér-közelítésben

A két rácshely kölcsönhatását leíró (Subscript[H, i, i + 1])Hamilton operátort a φ22 két-rácshely állapotok bázisán kellene tekinteni. Ezekre az elemi operátorokat szintén ezen a bázison kell definiálni. Ezt azonban most nem teszem meg, mert ez a tárgyalás nem szükséges az átlagtér-tárgyaláshoz. Ez utóbbi lényege éppen az, hogy a bonyolult, sok-rácshely problémát egy rácshely Hamilton operátorrá redukálja. Ennek az ára az, hogy az ilyen  (Subscript[H, MF])operátor új paramétereket tartalmaz, amelyeknek értékét önkonzisztencia-feltételekbõl kell meghatározni.
A két rácshely kölcsönhatását leíró Hamilton operátor legyen most az alábbi alakú:

h21 = {{b1, {sz}, {sz}}, {b2, {τy}, {τy}}};

Ezt a vektoros alakot fogom átalakítani átlagtér-Hamilton-operatorrá, amely csak egy-rácshely operátorokból áll: (Subscript[H, i, i + 1])->(Subscript[H, i]) + (Subscript[H, i + 1]). A fenti, vektoros alakban felírt operátor például (a b1, b2... együtthatók különbözõ csatolási paraméterek, amelyekre kesõbb feltételezéseket fogunk tenni): (Subscript[h2, i, i + 1]) = b1 (τx[i] ** 
                  sy[i + 1] ** ηy[i + 1] + τx[i + 1] ** 
                  sy[i] ** ηy[i] - τz[i] ** 
                  sx[i + 1] ** ηy[i + 1] - τz[i + 1] ** 
                  sx[i] ** ηy[i]) +b2 τx[
                i] ** τx[i + 1] + b3(
            sy[i] ** ηy[i] ** sy[i + 1] ** ηy[i + 1] + 
              sx[i] ** ηy[i] ** 
                sx[i + 1] ** ηy[i + 1])

Ezt a meanfield függvénnyel szeretném átalakítani (Subscript[H, i]) + (Subscript[H, i + 1])alakba, mégpedig úgy, hogy ahol az i. es az i+1. rácshelyen vett operátorok szorzata szerepel, ott egy-egy operátor esetében feltételezem, hogy a várható értéktõl való eltérése kicsi, és az ebben az eltérésben másodrendû tagokat elhanyagolom. A várható értéket a v függvény jelzi, ahol pl. v[{τx}] a τx operátor várható értéke a rendszer alapállapotában, (amit önkonzisztensen meg akarunk határozni), és pl. v[{τx,sy}] a  τx**sx operátor várható értéke, amely τx várh. értékétõl független.

ClearAll[kifbolop, mfkifejt, meanfield, aa, 
        bb];
kifbolop[kif_] := 
        Module[{bb = 1}, 
          For[k = 1, k <= Length[kif], 
            bb = bb ** kif[[k]][1]; k++]; 
          bb];
mfkifejt[aa_] := 
        2  aa[[1]]*
          If[SameQ[aa[[2]], 
              aa[[3]]], (-(((v[aa[[2]]]))^2)/2 + 
                kifbolop[aa[[2]]]*
                  v[aa[[2]]]), (-v[aa[[2]]]*v[aa[[3]]] + 
                v[aa[[2]]]*kifbolop[aa[[3]]] + 
                v[aa[[3]]]*
                  kifbolop[aa[[2]]])];
meanfield[h_] := Plus @@ Map[mfkifejt, h, {1}]

A mfkifejt és a kifbolop függvények mûködése egyszerûbb operátorokon demonstrálva:

aa = {b1, {τx}, {τx}}

mfkifejt[aa] // Simplify

kifbolop[{sy}]

Az átlagtér-megoldáshoz (amely a következõ lépésben egy-rácshely problémává alakul át) a h21 két-rácshely operátor átlagtér-kifejtett alakjához hozzáadjuk még a  λ*2τy[1]**sz[1] tagot ("spin-pálya kölcsönhatás") és két olyan tagot, amely a szabadenergia késõbbi deriválásához szükséges, és amelyek effektív mágneses illetve pszeudo-mágneses terek hatását írják le.

h2mf = (nce[meanfield[h21]] + 2  λ τy[1] ** sz[1] - 
            hz*sz[1] - hτy*τy[1]) // Simplify

Az átlagtér-Hamilton-op. mátrixa az egy-rácshely bázison, és ennek sajátértékei:

h2mat = 
        Table[nce[φa[[i]] ** h2mf ** φ[[j]]], {i, 
            8}, {j, 8}];

h2mat // MatrixForm

sert = Eigenvalues[h2mat] // Simplify

Eigensystem[h2mat] // Simplify

Termodinamika, a várható értékek kiszámítása

Termodinamikai mennyiségek definiálása a hõmérséklet függvényében

Az állapotösszeg (a benne szereplõ beta (=1/T) fordítottan arányos a hõmérseklettel, mint paraméterrel):

allapotosszeg = Plus @@ (Exp[-beta*sert])

Szabad energia számítása

fszaben = -Log[allapotosszeg]/beta;

Most feltételezéseket teszünk a csatolási együtthatókra, pl:

b1 = -13; b2 = -8; λ = 
      4;

A v várható értékeket szeretnénk meghatározni önkonzisztensen, ezzel megkapnánk a keresett átlagtér megoldást, mint a hõmérséklet(beta) függvényét. Az önkonzisztencia-feltételeket az adja, hogy az fszaben szabad energia bizonyos deriváltjai egyenlõk ezekkel a paraméterekkel (jelen esetben a hz "mágneses" tér szerinti derivált -1-szerese egyenlõ sz várható értékével, a hτy pályatérbeli "pszeudo-márgneses tér" szerinti derivált pedig τy várható értékét kell, hogy adja). Ezen egyenletek megoldása sokfajta megengedett várható érték esetén bonyolult lehet, de végsõ célunk ezen várható értékek meghatározása a hõmérséklet (1/beta) függvényében.

Az alábbi két kifejezés (vszfbol és vτyfbol) tehát v[{sz}]-vel ill. v[{τy}]-nal kell, hogy egyenlõ legyen.

ClearAll[beta]

ClearAll[vszfbol];
vszfbol[bet_, hzz_, htt_, vs_, vt_]:=-D[fszaben,  hz]/.{hτy -> htt,v[{sz}] -> vs,v[{τy}] -> vt,beta -> bet,hz -> hzz}//Simplify

ClearAll[vτyfbol];
vτyfbol[bet_, hzz_, htt_, vs_, vt_]:=-D[fszaben,  hτy]/.{hτy -> htt,v[{sz}] -> vs,v[{τy}] -> vt,beta -> bet,hz -> hzz}//Simplify

vτyfbol[bet, hzz, htt, vs, vt]
vszfbol[bet, hzz, htt, vs, vt]

Szuszceptibilitás, fajhõ: második deriváltak, definíció szerint:

ClearAll[χsz, χτy, 
        fajho];
χsz[bet_, hzz_, htt_, vs_, 
        vt_] := D[vszfbol[bet, hzz, htt, vs, vt], hzz] // 
        Simplify
χτy[bet_, hzz_, htt_, 
        vs_, vt_] := 
      D[vτyfbol[bet, hzz, htt, vs, vt], htt] // 
        Simplify
fajho[bet_, hzz_, htt_, vs_, vt_]:=-D[fszaben, {beta, 2}]/bet/.{hτy -> htt,v[{sz}] -> vs,v[{τy}] -> vt,beta -> bet,hz -> hzz}//Simplify

fajho[bet, hzz, htt, vs, vt]

Plot[vτyfbol[10, hz, 0.1, 0.2], {hz, 0.00001, 10}, 
      PlotRange -> All]

Az önkonzisztencia-egyenletek megoldása

hτy-ra ezután már nem lesz szükség, ezért 0-vá tesszük. 
Az egyszeruség kedvéért tegyük fel, hogy a külso mágneses tér is 0.

ClearAll[vszfbolh0, vτyfbolh0];
vszfbolh0[beta_, vs_, vt_] := 
        vszfbol[beta, hz, ht, vs, vt] /. {hz -> 0, 
            ht -> 0};
vτyfbolh0[beta_, vs_, vt_] := 
        vτyfbol[beta, hz, ht, vs, vt] /. {hz -> 0, 
            ht -> 0};

vτyfbolh0[beta, vs, vt]

Az explicit megoldás sajnos túl hosszú idõt venne igénybe,...

vszfbolh0[10, 0.1, 0.1] // N

... ezért numerikusan oldjuk meg a problémát. Ez azt jelenti, hogy különbözõ (1/)hõmérséklet értékeknél a vs és a vt értékeket ((Subscript[s, z]) illetve (Subscript[τ, y]) várható értékét) numerikusan határozzuk meg. Megkeressük a megoldást egy konkrét hõmérsékleten:

FindRoot[{vszfbolh0[1, vs, vt] == vs, 
        vτyfbolh0[1, vs, vt] == vt}, {vs, 0.4}, {vt, 0.4}]

FindRoot[{vszfbolh0[1, vs, vt] == vs, 
        vτyfbolh0[1, vs, vt] == vt}, {vs, 0.01}, {vt, 0.01}]

Mint látjuk, több megoldás is van, ezek közül az lesz a stabil, amelyiknek a hz ill. hτy szerinti deriváltja pozitív, ez további vizsgálatokat igényelne.

Keressük a megoldást a hõmérséklet (T) függvényében, tegyük fel, hogy 1/T az alábbi értékeket veheti fel:

inverzhomers = Range[0.01, 10. , 0.05];

mo = FindRoot[{vszfbolh0[#, vs, vt] == vs, 
              vτyfbolh0[#, vs, vt] == vt}, {vs, 0.3}, {vt, 0.3}] & /@ 
        inverzhomers

A kapott megoldás ábrázolása és kiértékelése

Ábrázoljuk a megoldást az (inverz) hõmérséklet függvényében

vsmo = Transpose[
            Partition[Flatten[Apply[List, mo, {2}]], 
              4]][[2]];
vtmo = Transpose[
            Partition[Flatten[Apply[List, mo, {2}]], 
              4]][[4]];

ListPlot[Abs[vsmo], PlotStyle -> {PointSize[0.02], Hue[ .6]}, 
      PlotRange -> All, 
      AxesLabel -> {"1/T", "v[{sz}]"}]
ListPlot[Abs[vtmo], PlotStyle -> {PointSize[0.02], Hue[ .6]}, 
      PlotRange -> All, 
      AxesLabel -> {"1/T", "v[{τy}]"}]

<< Graphics`MultipleListPlot`

MultipleListPlot[vsmo, vtmo, 
      SymbolStyle -> {GrayLevel[0.6], GrayLevel[0]}, 
      AxesLabel -> {"1/T", "(v[{sz}],v[{τy}])"}, 
      PlotLegend -> {"v[{sz}]", "v[{τy}]"}]

MultipleListPlot[vsmo, vtmo, 
      SymbolStyle -> {GrayLevel[0.6], GrayLevel[0]}, 
      PlotRange -> All, 
      AxesLabel -> {"1/T", "(v[{sz}],v[{τy}])"}, 
      PlotLegend -> {"v[{sz}]", "v[{τy}]"}]

homers = 1/inverzhomers;

vsmonew = 
        Transpose[{homers, Abs[vsmo]}];
vtmonew = Transpose[{homers, Abs[vtmo]}];

MultipleListPlot[vsmonew, vtmonew, 
      SymbolStyle -> {GrayLevel[0.6], GrayLevel[0]}, 
      PlotRange -> {{0, 20}, {0, 0.55}}, 
      AxesLabel -> {"T", "(v[{sz}],v[{τy}])"}, 
      PlotLegend -> {"v[{sz}]", "v[{τy}]"}]

MultipleListPlot[vsmonew, vtmonew, 
      SymbolStyle -> {GrayLevel[0.6], GrayLevel[0]}, 
      AxesLabel -> {"T", "(v[{sz}],v[{τy}])"}, 
      PlotLegend -> {"v[{sz}]", "v[{τy}]"}]

Látjuk, hogy a várható értékek egy kicsit zavarosan viselkednek alacsony hõmérsékleten, de magas T esetén 0-hoz tart az értékük, ahogyan várjuk (a hõmozgás lerombolja a rendezõdést). A "zavaros viselkedés" oka az lehet, hogy nem vizsgáltuk meg a kapott állapotok stabilitását, azaz, hogy a többféle megoldás közül melyiket kell választani. Ezt a hz ill. hτy szerinti újbóli deriválással kell meghatározni.

A megoldás stabilitásának vizsgálata. 
A stabil (szabadenergia minimumát biztosító) megoldás esetében a szuszceptibilitások és a fajhõ mindegyike pozitív kell, hogy legyen.

mofajho = 
      MapThread[
          fajho[#1, hzz, htt, #2, #3] &,  {inverzhomers, vsmo, 
            vtmo}] /. {hzz -> 0, htt -> 0}

moχsz = 
      MapThread[χsz[#1, hzz, htt, #2, #3] &,  {inverzhomers, vsmo, 
            vtmo}] /. {hzz -> 0, 
          htt -> 0}
moχτy = 
      MapThread[χτy[#1, hzz, htt, #2, #3] &,  {inverzhomers, vsmo, 
            vtmo}] /. {hzz -> 0, htt -> 0}

MultipleListPlot[moχsz, moχτy, mofajho, 
      SymbolStyle -> {GrayLevel[0.6], GrayLevel[0], GrayLevel[0.8]}, 
      AxesLabel -> {"1/T", "moχsz, moχτy, \
mofajho"}, 
      PlotLegend -> {"moχsz", "moχτy", "mofajho\
"}]

Látható, hogy jelen esetben mindhárom második derivált pozitív, ezzel a stabilitás vizsgálatát lezártam.

MultipleListPlot[moχsz, moχτy, 
      SymbolStyle -> {GrayLevel[0.6], GrayLevel[0]}, 
      AxesLabel -> {"1/T", "moχsz, moχτy"}, 
      PlotLegend -> {"moχsz", "moχτy"}, 
      PlotRange -> {{0, 20}, {0, 0.2}}]

A szuszceptibilitásokban látható csúcsok ill. ugrások elvileg fázisátalakulást jeleznek az anyagban, azon a hõmérsékleten, ahol az ugrás található. Ezeket közelebbrõl, nagyobb felbontásban vizsgálva, többet is lehet késõbb mondani róluk.

Irodalom

Apagyi B.: Kvantummechanika, Mûegyetemi Kiadó, 1999

Sólyom J.: A modern szilárdtestfizika alapjai I., II., III., ELTE Eötvös Kiadó, 2003