Skip to main content

Modelling quantum spin-system with orbital degeneracy

· 6 min read
Kirill Vasin

This notebook delves into quantum mechanics with a focus on symbolic quantum operators and spin-Hamiltonians. It begins with setting up orthogonality rules and defining the linear properties of basis functions.

Next, we investigate the spin-Hamiltonian for a cubic system with doubly degenerated orbital states, particularly focusing on Fe2+Fe^{2+} ions in cubic spinels with 5Eg^5E_g states. The Hamiltonian for a paramagnetic state is presented and implemented in code, followed by defining basis vectors and functions to find eigenvalues and eigenvectors, and calculating thermal averages of quantum operators.

Let us start from the ground level - how to represent a state

Download original notebook 
(*BB[*)(Ket[(*FB[*)((1)(*,*)/(*,*)(2))(*]FB*)])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)

Just copy and paste the expression above to have prettier representation of your wavefunction.

Let's work with wavefunction in the basis of S=2S=2 and a doubly-degenerated orbital state L^=1/2\hat{L}=1/2

(*BB[*)(Ket[Ms, Ml])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)

Setting up rules​

Firstly, define the ortogonality of our basis functions

Ket /: (*BB[*)(Bra[Ms1_, Ml1_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*).(*BB[*)(Ket[Ms2_, Ml2_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*) := KroneckerDelta[Ms1, Ms2] KroneckerDelta[Ml1, Ml2]

Ket /: Conjugate[(*BB[*)(Ket[Ms2_, Ml2_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)] := (*BB[*)(Bra[Ms2, Ml2])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*)

Linear properties must also be defined as well

Unprotect[Dot];
Unprotect[Conjugate];

Conjugate[a_ b_] := Conjugate[a] Conjugate[b]
Conjugate[a_ + b_] := Conjugate[a] + Conjugate[b]
Dot[a_, b_?NumericQ] := b a
Dot[b_?NumericQ, a_] := b a
Dot[a_ + b_, c_] := a.c + b.c 
Dot[c_, a_ + b_] := c.a + c.b  
Dot[a_, b_?NumericQ c_] := b a.c
Dot[b_?NumericQ c_, a_] := b c.a

Protect[Dot];
Protect[Conjugate];

Symbolic quantum operators​

A well known ladder-like operators (see more here) as well as pauli-like operators UθU_\theta, UεU_\varepsilon for an orbital doublet

Ket /: Sp.(*BB[*)(Ket[Ms_, Ml_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*) := (*SqB[*)Sqrt[2(2+1) - (Ms+1)Ms](*]SqB*) (*BB[*)(Ket[Ms+1, Ml])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)
Ket /: Sm.(*BB[*)(Ket[Ms_, Ml_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*) := (*SqB[*)Sqrt[2(2+1) - (Ms-1)Ms](*]SqB*) (*BB[*)(Ket[Ms-1, Ml])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)
Ket /: Sz.(*BB[*)(Ket[Ms_, Ml_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*) := Ms (*BB[*)(Ket[Ms, Ml])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)

Uθ = Sum[(*BB[*)(Ket[Ms,1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*).(*BB[*)(Bra[Ms, 1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*) - (*BB[*)(Ket[Ms,-1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*).(*BB[*)(Bra[Ms, -1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*), {Ms,-2,2}];
Uϵ = Sum[(*BB[*)(Ket[Ms,1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*).(*BB[*)(Bra[Ms, -1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*) + (*BB[*)(Ket[Ms,-1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*).(*BB[*)(Bra[Ms, 1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*), {Ms,-2,2}];

Sx = (*FB[*)((1)(*,*)/(*,*)(2))(*]FB*)(Sp + Sm);
Sy =  (*FB[*)((1)(*,*)/(*,*)(2I ))(*]FB*)(Sp - Sm);

Testing​

check a single operator

Sm.(*BB[*)(Ket[0,1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)
((*SqB[*)Sqrt[6](*]SqB*)) ((*BB[*)(Ket[-1,1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*))

We can also build the matrix of a operator in the basis of spin-only states

Table[

  (*BB[*)(Bra[Ms1, 1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*).(Sz).(*BB[*)(Ket[Ms2, 1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)

, {Ms1, -2,2}, {Ms2, -2,2}] // MatrixForm 
((*GB[*){{-2(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)-1(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)1(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)2}}(*]GB*))

Now check the modulus of an spin-momentum

Table[

  (*BB[*)(Bra[Ms1, 1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*).(Sz.Sz + Sx.Sx + Sy.Sy).(*BB[*)(Ket[Ms2, 1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)

, {Ms1, -2,2}, {Ms2, -2,2}] // MatrixForm 
((*GB[*){{6(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)6(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)0(*|*),(*|*)6(*|*),(*|*)0(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)6(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)6}}(*]GB*))

This is correct! As it should be S(S+1)=6S(S + 1) = 6

Spin-Hamiltonian​

Let us have a look at the cubic system with a doubly degenerated orbital states. For example, it can be effectively Fe2+Fe^{2+} ions in a cubic spinels with 5Eg^5E_g states as a ground multiplet. Let us have a look at the Hamiltonian for a paramagnetic state

H=μg SH0−d[Uθ(3Sz2−S(S+1))+3Uϵ(Sx2−Sy2)]H = \mu g~ \mathbf{S H_0} - d \Big[ U_{\theta} (3 Sz^2 - S(S+1)) + \sqrt{3} U_{\epsilon}(Sx^2 - Sy^2) \Big]

where dd accounts for an effective spin-orbit coupling and then dictates cubic anistoropy for the spin

Implementation​

Now we move to the actual code

μ = 0.466;
g = 2.;

Effective hamiltonian in our notations

H[Hx_, Hy_, Hz_, d_] := μ g (Hx Sx + Hy Sy + Hz Sz ) - d (3 Sz.Sz.Uθ -  2(2+1) Uθ + (*SqB[*)Sqrt[3](*]SqB*) Sx.Sx.Uϵ - (*SqB[*)Sqrt[3](*]SqB*) Sy.Sy.Uϵ)

Now define an array of basis vectors and function to find eigenvalues and eigenvectors

basis = Table[(*BB[*)(Ket[Ms2, Ml])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*), {Ms2,-2,2}, {Ml,-1,1,2}] // Flatten; 

system[params__] := With[{ h = Simplify[H[params]]}, With[{
  r = Table[

    Conjugate[i].(h.j)

  , {i, basis}, {j, basis}]  // Eigensystem // Transpose 
},

  (* truncate accuracy *)
  r // Chop
] ]

Now we need a helper function to find a thermal average of any quantum operator

(* for caching already calculated values *)

ClearAll[cached];
cached[expr_] := cached[expr] = expr;
SetAttributes[cached, HoldFirst]


average[sys_, T_, operator_] := With[{
  operatorMatrix = cached[Table[Conjugate[i].operator.j, {i, basis}, {j, basis}]],
  
  stat = Total[Exp[-(*FB[*)((#[[1]])(*,*)/(*,*)(0.7 T))(*]FB*)] &/@ sys]
},
  Map[
    ((Conjugate[{#[[2]]}].operatorMatrix.Transpose[{#[[2]]}]) // First // First) (*FB[*)((1)(*,*)/(*,*)(stat))(*]FB*) Exp[-(*FB[*)((#[[1]])(*,*)/(*,*)(0.7 T))(*]FB*)] &
    , sys] // Total  
]

Usually for such systems a parameter d≃3−4 cm−1d\simeq 3-4~cm^{-1}.

Let's try to see a magnetic anisotropy in an averaged magnetic moment of a single ion site at T=10 KT=10~K in the magnetic field around 10 T10~T

Table[
  With[{
    (*BB[*)(*apply a magnetic field to the given direction*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
    sys = system[10 Sin[x], 10 Cos[x], 0, 3.]
  },
    (*BB[*)(*project magnetic moment to the direction of field*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
    (2.0 average[sys, 10, #] &/@ {Sx, Sy}).({Sin[x],  Cos[x]})
   
  ] {Sin[x],  Cos[x]}
, {x, 0, 2Pi, 0.1}] // Re;

ListLinePlot[%, 
  PlotRange->Full, 
  Epilog->{
    Black // Lighter, Arrow[{{0, 0}, {2.5501, 0.0}}],
    Arrow[{{0,0}, {0, 2.52099}}]
}, AxesLabel->{"\\mu x", "\\mu y"}, ImageSize->500, AspectRatio->1]
(*VB[*)(FrontEndRef["c7710bc3-f3fa-4da9-b610-78b5b046c42c"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJ5ubGxokJRvrphmnJeqapCRa6iaZGRromlskmSYZmJglmxglAwCKFxXG"*)(*]VB*)

It is clear, that a cubic anistropy is quite visible here 💫

Actually, this is not a typical spin system with a fixed magnetic anisotropy. Here it comes from the orbital degeneracy, which is lifted in magnetic fields. However, it also means, that the field can change anisotropy. The less field we apply, the less spin prjection we have, the less splitting between the orbital states we have, and therefore the less anisotropy we observe. Try to evaluate the cell above using different fields

Ferromagnetic ordering​

One can go further and model an ordered state using mean field approximation, which is suitable for ferromagnets

H=μg SH0−d[Uθ(3Sz2−S(S+1))+3Uϵ(Sx2−Sy2)]+JSMH = \mu g~ \mathbf{S H_0} - d \Big[ U_{\theta} (3 Sz^2 - S(S+1)) + \sqrt{3} U_{\epsilon}(Sx^2 - Sy^2) \Big] + J \mathbf{S M}

where MM is an average of S^\hat{S} (i.e. just a classical vector).

Implementation​

H[Hx_, Hy_, Hz_, d_, J_, Mx_, My_, Mz_] := μ g (Hx Sx + Hy Sy + Hz Sz ) - d (3 Sz.Sz.Uθ -  2(2+1) Uθ + (*SqB[*)Sqrt[3](*]SqB*) Sx.Sx.Uϵ - (*SqB[*)Sqrt[3](*]SqB*) Sy.Sy.Uϵ) + J (Mx Sx + My Sy + Mz Sz)

Now it becomes more computationally expensive to do. It is better monitor the progress, for this reason we need a helper function


Progress bar​

progressBar.wl

Module[{progressBar},
  progressBar[max_, OptionsPattern[] ] := With[{cell = OptionValue["Cell"]}, LeakyModule[{progress = 0., indicator = 0., handler, timer = AbsoluteTime[]},
    With[{
        cell = CellPrint[ToString[
            Graphics[{
              LightBlue, 
              Rectangle[{-0.004,-0.12}, {1.004,1.12}], Green, 
              Rectangle[{0,0}, {Offload[indicator],1}]
            }, PlotRange->{{0,1}, {0,1}}, ImagePadding->0, ImageSize->{300, 20}, Controls->False]
        , StandardForm], "After"->cell]
    },
        handler[n_] := With[{diff = AbsoluteTime[] - timer},
            progress = n;
            If[diff > 0.03,
                indicator = progress / max // N;
                timer = AbsoluteTime[];
            ];

            If[Abs[progress - max] < 0.1, Delete[cell] ];
        ];

        handler
    ]
  ] ];

    Options[progressBar] = {"Cell":>EvaluationCell[]};

    progressBar
]

we store it to a file, so we can reuse it later

Continue with our plotting

Since this is a mean approximation method, we need to substitute an averaged value of a spin back to the Hamiltonian and wait until it converges. The memory of our system is stored exactly in this mean field vector

meanField = {0.1, -0.8, -0.1};
sys;

progressBarGenerator = Get["progressBar.wl"];

With[{
  progress = progressBarGenerator[2Pi]
},

Table[
  With[{
    
  },
  
    Do[
    
      sys = system @@ Join[{10. Sin[x], 10. Cos[x], 0, 3., -100.}, meanField];
      meanField = average[sys, 10, #] &/@ {Sx, Sy, Sz}; 
      
    , {i, 15}];

    progress[x];
    
    
    ( meanField[[{1,2}]] ) . {Sin[x], Cos[x]}
   
  ] {Sin[x], Cos[x]}
  
, {x, 0, 2Pi, 0.05}]

];

ListLinePlot[%, 
  PlotRange->Full, 
  Epilog->{
    Black // Lighter, Arrow[{{0, 0}, {2.3501, 0.0}}],
    Arrow[{{0,0}, {0, 2.32099}}]
}, AxesLabel->{"\\mu x", "\\mu y"}, ImageSize->500, AspectRatio->1]
(*VB[*)(FrontEndRef["8d6b2390-be0d-45a5-b9d7-099c0e12ef55"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKW6SYJRkZWxroJqUapOiamCaa6iZZppjrGlhaJhukGhqlppmaAgCErRWn"*)(*]VB*)

From this plot one can see, that our magnetic moment sticks to C4C_4 directions and when the rotating magnetic field overcomes the anisotropy magnetic moment switches to another easy-axis or direction.

A calculation itself is quite slow, due to algebra we used for quantum computations. In general, we do not need all operators, once we constructed the Hamiltonian in the basis we have. We can deal only with matrixes and simple numerics.

.md
### Converting operators to matrixes
We build a function, that approximate any arbitary function up to 2nd order for all variables in arguments

Converting operators to matrixes​

We build a function, that approximate any arbitary function up to 2nd order for all variables in arguments

quadInterpolator.wl

BeginPackage["Misc`QuadarticInterpolation`"]
Begin["`Private`"]

quadraticInterpolation[f_, vars__, OptionsPattern[]] := Module[{
  variables =  List[vars],
  rules,
  constant,
  quadratic,
  final,
  cached,
  mdcache,
  calculated = <||>,
  linear,
  handler, count = 0
},

  If[OptionValue["Monitor"] =!= Null,
    handler = OptionValue["Monitor"][ (Length[variables] Length[variables] - Length[variables]) / 2  + 2 Length[variables]];
  ];

  cached[expr_, {i_,j_}] := mdcache[expr, Sort[{i,j}]];
  SetAttributes[cached, HoldFirst];

  mdcache[expr_, {i_,j_}] := If[KeyExistsQ[calculated, {i,j}],
    calculated[{i,j}]
  ,
    calculated[{i,j}] = expr
  ];

  SetAttributes[mdcache, HoldFirst];

  rules = (#[[1]] -> #[[2]]) &/@ variables;
  
  constant = Hold[f] /. rules // ReleaseHold;

  linear = Table[
    With[{
      r = {i[[1]] -> 1.}
    },

      count++;
      handler[count];

      (Hold[f] /. r /. rules // ReleaseHold) - constant
    ]
    
  , {i, variables}];

  quadratic = Table[
    cached[With[{
      f11 = With[{
        r = {i[[1]] -> 1., j[[1]] -> 1.}
      },
        (Hold[f] /. r /. rules // ReleaseHold)
      ],

      fm1m1 = With[{
        r = {i[[1]] -> -1., j[[1]] -> -1.}
      },
        (Hold[f] /. r /. rules // ReleaseHold)
      ],

      f1m1 = With[{
        r = {i[[1]] -> 1., j[[1]] -> -1.}
      },
        (Hold[f] /. r /. rules // ReleaseHold)
      ],

      fm11 = With[{
        r = {i[[1]] -> -1., j[[1]] -> 1.}
      },
        (Hold[f] /. r /. rules // ReleaseHold)
      ]      
      
    },
      count++;
      handler[count];

      Print["Checking... ", i, j];

      (*FB[*)((f11 - f1m1 - fm11 + fm1m1)(*,*)/(*,*)(4.))(*]FB*)
    ], {i,j}]
    
  , {i, variables}, {j, variables}];

  final = Sum[
    With[{l = linear[[i]]},
      Slot[i] Hold[l]
    ]
  , {i, Length[linear]}] + Sum[
    With[{
      q = (*FB[*)((1)(*,*)/(*,*)(2))(*]FB*)quadratic[[i, j]]
    },
      Slot[i]Slot[j] Hold[q]
    ]
  , {i, Length[linear]}, {j, Length[linear]}] + With[{c = constant}, Hold[c] ];

  With[{
    r = final
  },
    (r &) /. {Hold -> Identity}
  ]
]

SetAttributes[quadraticInterpolation, HoldFirst]

Options[quadraticInterpolation] = {"Monitor"->Null}

End[]
EndPackage[]

Misc`QuadarticInterpolation`Private`quadraticInterpolation

Now lets import it

quadInterpolation = Get["quadInterpolator.wl"];
progressBarGenerator = Get["progressBar.wl"];

Now we apply it to our hamiltonian calculated in the basis with a given parameters

int = quadInterpolation[

  Table[Conjugate[fin].H[a1,a2,a3, a4, a5, a6,a7,a8].ini, {fin, basis}, {ini, basis}]
  
, {a1,0}, {a2,0}, {a3,0}, {a4,0}, {a5,0}, {a6,0}, {a7,0}, {a8,0}, "Monitor"->progressBarGenerator];

Now one can check that it speeds up things drastically

(* old *)
Table[Conjugate[fin].H[0,1.,1., 1., 1.,  1.,1.,1.].ini, {fin, basis}, {ini, basis}] // RepeatedTiming // First 
(* converted *)
int[0,1.,1., 1., 1.,  1.,1.,1.] // RepeatedTiming // First 
1.051729`
0.0007063037109375`

We can also save this function for later, so that we do not need to derive it again

int // Iconize 
(*VB[*)(Get[FileNameJoin[{".iconized", "iconized-6f9a.wl"}]])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRAeH5DwTM7Py6xKdcvMSXXKr8i0EGdmAADrVArH"*)(*]VB*)

Or better to internal Notebook storage, so that we do not need to worry about folders structure

NotebookStore["Int"] = int
"Int"

Modelling​

using faster approach

progressBarGenerator = Get["progressBar.wl"];
int = NotebookStore["Int"];

meanField = {0.1, -0.8, -0.1};
sys;

With[{
  progress = progressBarGenerator[15]
},
Table[
  With[{
    
  },
  
    Do[
      sys = int @@ Join[{h Sin[x], h Cos[x], 0, 3., -100.}, meanField] // Eigensystem // Transpose;
      meanField = average[sys, 10, #] &/@ {Sx, Sy, Sz}; 
    , {i, 25}];

    progress[h];

    
    
    ( meanField[[{1,2}]]).({Sin[x], Cos[x]})
   
  ] {Sin[x], Cos[x]}
, {h, {1,5,10,15}}, {x, 0, 2Pi, 0.05}]
];

ListLinePlot[%, 
  PlotRange->Full, 
  Epilog->{
    Black // Lighter, Arrow[{{0, 0}, {2.5501, 0.0}}],
    Arrow[{{0,0}, {0, 2.52099}}]
}, AxesLabel->{"\\mu x", "\\mu y"}, ImageSize->500, AspectRatio->1, PlotLegends->{"1T", "5T", "10T", "15T"}]
(*VB[*)(FrontEndRef["17a8d4a4-9c68-4b79-a0ce-37fa69f6fee7"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKG5onWqSYJJroWiabWeiaJJlb6iYaJKfqGpunJZpZppmlpaaaAwCGQBYd"*)(*]VB*)

Here we have magnetic moment curves plotted for different magnitudes of applied field

(*GB[*){{(*VB[*)(RGBColor[0.368417, 0.506779, 0.709798])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRcVTJv0VGXKdfsiw2kvp3eYP7AvOt4yw3vVtmf2AKMBH1E="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.880722, 0.611041, 0.142051])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRf5HIl+fl/vjX3Rv+oPt5b2PrYv2jxrnfouvUP2ALwBIFo="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.560181, 0.691569, 0.194885])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRcxPIvawvD2oX1R6dq9/qGKz+yLJCfsvff34wl7AJc8HzE="*)(*]VB*)(*|*),(*|*)(*VB[*)(RGBColor[0.922526, 0.385626, 0.209179])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeGJAIcndyzs/JLwouTyxJzghJzS3ISSxJTWMGyXMgyRdd2frNLrT1rX3RHf35ayTW3bAv0nBnXJRw7JQ9AKbfHmU="*)(*]VB*)}(*||*),(*||*){1(*|*),(*|*)5(*|*),(*|*)10(*|*),(*|*)15}}(*]GB*)

Why do we have circles at the bottom?

Well, it looks that the magnetic field cannot overcome the anisotropy and the position of magentic moment remains constant, while projection goes from 0 up to the maximum. We can prove its shape by plotting this function

PolarPlot[(*BB[*)(*our constant magnetic moment*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) {0,-1}.{ (*BB[*)(*projection*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*) Cos[x], Sin[x]}, {x,0,2Pi}, PlotRange->{{-1,1},{-1,1}}, Epilog->{SVGAttribute[Circle[{0,0}, 1], "stroke-dasharray"->"3"]}] 
(*VB[*)(FrontEndRef["80ee3f9f-f7c6-4c81-9107-63b95fea7fa0"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWxikphqnWabpppknm+maJFsY6loaGpjrmhknWZqmpSaapyUaAACJHBXn"*)(*]VB*)

Real-time modelling​

We can go further and recalculate spins based on user's input. For this case Offload keyword is used together with some GUI building blocks

Download this notebook to try it

int = NotebookStore["Int"];

LeakyModule[{
  spin = {0.1, -0.8, -0.11},
  defaults = {0.,0.,0., 3., -100.},
  temperature = 10.,
  field = {0.,0.},
  controlPanel,
  frame = CreateUUID[]
},

  controlPanel = InputGroup[{
    InputRange[0, 5., 1, 3, "Label"->"Spin-orbit"],
    InputRange[-100, 0., 10, -100, "Label"->"J"],
    InputRange[1, 150., 10, 10, "Label"->"Temperature"]
  }];

  EventHandler[frame, Function[
    With[{
      sys = int @@ Join[defaults, spin] // Eigensystem // Transpose
    },
      spin = average[sys, temperature, #] &/@ {Sx, Sy, Sz} // Chop // Re;
    ]  
  ]];

  EventHandler[controlPanel, Function[data,
    defaults[[4;;5]] = data[[1;;2]];
    temperature = data[[3]];
  ]];

  {
    Graphics[{
      White, EventHandler[Rectangle[{-4,-4}, {4,4}], {"mousemove" -> Function[xy,
        defaults[[1;;2]] = - 10 xy;
        field = xy 0.8;
      ]}],
      Red, Arrow[{{0,0}, With[{s = spin}, {s[[1]], s[[2]]}] // Offload}],
      Blue, Arrow[{{0,0}, field // Offload}],
      AnimationFrameListener[spin // Offload, "Event" -> frame]
    }, PlotRange->1.5{{-2,2}, {-2,2}}, Controls->False, TransitionType->"Linear", TransitionDuration->30, ImagePadding->0],
    
    controlPanel
  } // Row  
]

Expected result is shown in this video