Spin expectation values

[thanasis]

I was wondering if there is a way to extract the spin expectation values of a multielectron system starting from a description of its Sys.S and Sys.ee.

I suppose that since ES carries out a full matrix calculation the information must be somewhere in there. Do I understand this correctly?

Thanks in advance!

[joscha_nehrkorn]

I have some good news: Yes, you can get any expectation value for any spin-operator based observable.
But, I also have bad news: You have to build it!

With

Code: Select all

H = sham(Sys,[0,0,0]);
[Vt,E]=eig(H);

you get eigen vectors in Vt and eigen values, i.e. energies, on the diagonal of E.

Code: Select all

Vt(:,1)

is your ground state.

Then you need to build operators with

Code: Select all

sop

.
S^2 ( I guess this is what you want to have, right?) for a 2-spin system would be:

Code: Select all

strxyz ={'x','y','z'};
s2t = zeros((2*Sys.S(1)+1)*(2*Sys.S(2)+1));
for n =3:-1:1
  for m=2:-1:1
    str = [strxyz{n}, num2str(m)];
    op(m,n,:,:)=sop(Sys,str);
  end
  opt(n,:,:) = squeeze(op(1,n,:,:)+op(2,n,:,:));
  s2t = s2t + squeeze(opt(n,:,:))*squeeze(opt(n,:,:));
end

<GS|S^2|GS> is than:

Code: Select all

Vt(:,1).'*s2t*Vt(:,1)

Just a warning: This will not give the spin! In cases with S a good quantum number you get S(S+1).

[thanasis]

Dear Joscha,

Thanks for this very didactic reply!

I suspected that I would need to do a bit of work myself, so the news are not really all that bad!

So, just to see if I understood everything correctly...
-s2t = zeros((2*Sys.S(1)+1)*(2*Sys.S(2)+1)); creates an empty square matrix at the dimensions of the Hilbert space of the problem
-The for n =3:-1:1 loop cycles through the x, y and z components and the for m=2:-1:1 subloop cycles through each of the spins
-Then with the sop function we calculate the Six, Siy and Siz operators for each spin

E.g., for S1 = S2 = 1/2, a ferromagnetic interaction gives a result of 2 (S = 1) and an antiferromagnetic a result of 0 (S = 0). For S1 = S2 = 1 I get 6 (S = 2) and 0 (S = 0), respectively, which make sense.

I am trying to understand what the squeeze matlab function does and how to use it for my purposes. in particular, I am trying to calculate the <Siz> (and subsequently <Six>, <Siy>) for each individual spin of a trinuclear cluster with S1 = S2 = S3 = 5/2. So, I suppose I should change your input to something like:

Code: Select all

strxyz ={'x','y','z'};
s2t = zeros((2*Sys.S(1)+1)*(2*Sys.S(2)+1)*(2*Sys.S(3)+1));
for n =3:-1:1
  for m=3:-1:1
    str = [strxyz{n}, num2str(m)];
    op(m,n,:,:)=sop(Sys,str);
  end
  opt(n,:,:) = squeeze(op(1,n,:,:)+op(2,n,:,:)+op(3,n,:,:));
  s2t = s2t + squeeze(opt(n,:,:))*squeeze(opt(n,:,:));
end

(although I am still not clear on the squeeze function).

Then how should I go about to calculate <Six>, <Siy> and <Siz> for each spin?

[thanasis]

As a side note, to get the actual S values, I used:

Code: Select all

SStot = Vt(:,1).'*s2t*Vt(:,1); % S(S+1)
allroots = roots([1 1 -SStot]); % Two roots
Stot = allroots(allroots>0); % S

which keeps only the positive root, since the 2nd degree polynomial x^2 + x - a = 0 (a > 0) always has one positive and one negative solution.

[joscha_nehrkorn]

Code: Select all

sop

do not calculate anything! It returns the corresponding operator. This is a very important difference.
Otherwise you understood the code correctly.

Code: Select all

squeeze

does remove singleton dimension (see Matlab help). Roughly speaking this means unnecessary dimension, as it is only one element. Here it is used to take care that only matrices are added and multiplied.
The operators for Six and so on are all contained in

Code: Select all

op

A small remark: In most cases EasySpin handle degeneracy in a clever way, but not always. In zero field for 3 half-integer spin you will always get degenerate states. Then also linear combination of the two are a solution, which might have unexpected expectation values of the corresponding operators. A small field (think about the direction) do help here a lot.

[thanasis]

Thanks Joscha,

I'm still trying to wrap my mind around the 4D matrices that sop produces, but until then I did a little hack. It isn't as elegant as your code, but I think it gives correct results:

Code: Select all

[S1x,S1y,S1z,S2x,S2y,S2z,S3x,S3y,S3z] = sop(Sys,'x1','y1','z1','x2','y2','z2','x3','y3','z3')
S1zvalue = Vt(:,1).'*S1z*Vt(:,1); % <S1z>
S2zvalue = Vt(:,1).'*S2z*Vt(:,1); % <S2z>
S3zvalue = Vt(:,1).'*S3z*Vt(:,1); % <S3z>

If I implemented this correctly, these are the spin expectation values <Siz> for each of my spins. Their sum is -1/2, which is in line with a ms = -1/2 ground state. I get the same result for many different Jij combinations and also for setting H = sham(Sys(i),[0,0,1]) (a small magnetic field along the z axis).

I hope I got this right. I will be back as soon as I have done some more tests!

[thanasis]

I am coming back to this question after doing some tests:
I managed to get things to work according to Joscha's advice. To sum up, the code I am using is:

Code: Select all

cm = 100*clight/1e6; % Conversion constant from cm-1 to MHz
Si = 1/2; % = 1/2, 3/2, 5/2
J12 = -20; J13 = -20; J23 = -25; % Isotropic exchange
dzz = 1.0; % Antisymmetric exchange z-component
Gx = 0; Gy = 0; Gz = dzz; % Antisymmetric tensor elements
GG = [0 Gz -Gy; -Gz 0 Gx; Gy -Gx 0]; % Antisymettric exchange tensor

Sys.S = [Si Si Si];
Sys.g = [gxy gz; gxy gz; gxy gz];
Sys.ee = cm*[-2*J12(i)*eye(3) - 2*GG; -2*J13(i)*eye(3) + 2*GG; -2*J23(i)*eye(3) - 2*GG];

H = sham(Sys(i),[0,0,0.1]); % Generate the Hamiltonian matrix of the system
[Vt,E]=eig(H); % Generate the right-hand vectors (Vt = |>) and the energy eigenvalues (E) of the Hamiltonian

strxyz ={'x','y','z'};
s2t = zeros((2*Sys(i).S(1)+1)*(2*Sys(i).S(2)+1)*(2*Sys(i).S(3)+1));
for l=3:-1:1 % Cycle trhough x (l = 1) / y (l = 2) / z (l = 3)
  for m=3:-1:1 % Cycle through spins 1/2/3
    str = [strxyz{l}, num2str(m)]; % Concatenate a string like x1, x2, x3, y1,...
    op(m,l,:,:)=sop(Sys(i),str);
  end
  opt(l,:,:) = squeeze(op(1,l,:,:)+op(2,l,:,:)+op(3,l,:,:));
  s2t = s2t + squeeze(opt(l,:,:))*squeeze(opt(l,:,:));
  end

[S1x,S1y,S1z,S2x,S2y,S2z,S3x,S3y,S3z] = sop(Sys,'x1','y1','z1','x2','y2','z2','x3','y3','z3');
S1zvalue = Vt(:,1).'*S1z*Vt(:,1);
S2zvalue = Vt(:,1).'*S2z*Vt(:,1);
S3zvalue = Vt(:,1).'*S3z*Vt(:,1);

As you may see, this is a triangle of half-integer spins with an antisymmetric exchange dzz parallel component.

In the special case where dzz = 0, the <Siz> values are all reals and their sum is -1/2, equal to the Ms value of the ground state. However, when I tried to introduce AE, I observed that the expectation values were becoming complex numbers.
Also, when I tried:

Code: Select all

S1xvalue = Vt(:,1).'*S1x*Vt(:,1);
S2xvalue = Vt(:,1).'*S2x*Vt(:,1);
S3xvalue = Vt(:,1).'*S3x*Vt(:,1);

S1yvalue = Vt(:,1).'*S1y*Vt(:,1);
S2yvalue = Vt(:,1).'*S2y*Vt(:,1);
S3yvalue = Vt(:,1).'*S3y*Vt(:,1);

I also got complex numbers even if dzz = 0.

Then I went to brush up my quantum mechanics and remembered that the <Six>, <Siy> and <Siz> are <GS|S2t|GS>, where <GS| is the complex conjugate.

So, I changed to:

Code: Select all

[S1x,S1y,S1z,S2x,S2y,S2z,S3x,S3y,S3z] = sop(Sys(i),'x1','y1','z1','x2','y2','z2','x3','y3','z3');
S1zvalue(i) = conj(Vt(:,1).')*S1z*Vt(:,1);
S2zvalue(i) = conj(Vt(:,1).')*S2z*Vt(:,1);
S3zvalue(i) = conj(Vt(:,1).')*S3z*Vt(:,1);

S1xvalue(i) = conj(Vt(:,1).')*S1x*Vt(:,1);
S2xvalue(i) = conj(Vt(:,1).')*S2x*Vt(:,1);
S3xvalue(i) = conj(Vt(:,1).')*S3x*Vt(:,1);

S1yvalue(i) = conj(Vt(:,1).')*S1y*Vt(:,1);
S2yvalue(i) = conj(Vt(:,1).')*S2y*Vt(:,1);
S3yvalue(i) = conj(Vt(:,1).')*S3y*Vt(:,1);

Now all expectation values are real numbers, as expected. Have I got it correct?

Also, I would like to note that when I apply a minute magn. field along the z-axis, <Six> and <Siy> are 12 to 13 orders of magnitude smaller than <Siz>. I was expecting that AE would cant the spins a little further away from the z-axis. What do you think?

[joscha_nehrkorn]

There was a mistake in my code.
The correct form is:

Code: Select all

Vt(:,1)'*S1z*Vt(:,1)

instead of

Code: Select all

Vt(:,1).'*S1z*Vt(:,1)

.
The latter is only the transpose, while the correct one is the complex conjugated transpose or Hermitian transpose.

Code: Select all

Vt(:,1)'

is identical to

Code: Select all

conj(Vt(:,1).')

.

Sorry for that!

There is quite a strong effect of the anisotropic interaction! Set isotropic interaction to zero and inspect the numbers.
In this case

Code: Select all

Vt(:,1)'*S1x*Vt(:,1)

and

Code: Select all

Vt(:,1)'*S1y*Vt(:,1)

are indeed numerical zeros and

Code: Select all

Vt(:,1)'*S1z*Vt(:,1)

is not 0.5! You see that you deal with quantum mechanical objects and not arrows, who point somewhere. It is not simple to identify the proper operator to make effects on the wave function more obvious.
Sz as the sum of the the three single spin operators is something that might make more sense to look at.

[thanasis]

Well put Joscha... a very didactic example.

I did "know" that spins are not classical arrows, but without a real model and calculations that hadn't really sunk in. To make matters worse, I have been reading quite a lot of PRL/PRB-level publications with nicely drawn arrows for spins "pointing" at various directions within the cluster. I was hoping to understand what the physical meaning of these directions is and I (naively?) thought that this could be solved by finding the <Sx/y/z> values and using them to draw vectors. Probably an oversimplification from my part... (I am more of an experimental chemist and not a theoretician after all).

What I did find by playing around with parameters is that <Sx> and <Sy> become significant for magnetic fields with transverse components. Axial magn. fields only produce a significant <Sz> value.

Anyway, I'll try to understand this better.

[Stefan Stoll]

It is absolutely legitimate to draw arrows representing the the expectation value of any vector spin operator, for example <Svec>_q = [<q|Sx|q>;<q|Sy|q>;<q|Sz|q>] for some eigenstate |q>, as long as it is explicitly and clearly stated what the arrow represent. These arrows can help rationalize behavior of spin systems. One just needs to be careful to always be aware what the arrows represent.