Bug for coupled spins in pepper near anticrossings?

[theresia_q]

When using the ‘pepper’ function, very strange spectra occur for coupled systems in the vicinity of anticrossings between the levels. As an example, I have simulated two such spectra, one for a doublet-doublet pair and one for a triplet-doublet pair:

ddp.png (42.82 KiB) Viewed 13333 times

tdp.png (54.55 KiB) Viewed 13333 times

It doesn't matter whether you use ‘initState’ (and the basis doesn't matter either) or a Boltzmann population. The code that creates the example spectra is as follows:

Code: Select all

%% doublet doublet pair
clear all

% Spin system
DDP.J = 9500; % Change the exchange interaction, a value was selected here so that anticrossings appear

DDP.S = [1/2 1/2];
DDP.g = [2.0020; 1.9999];
DDP.dip = 350;

DDP.lwpp = 2; 
%RP.initState = 'singlet'; % switch on and off spin polarisation

% Experimental parameters
Exp.Range = [310.777, 374.777];
Exp.mwFreq = 9.6;
Exp.Harmonic = 0;
Exp.Temperature = 80;


[field, sim] = pepper(DDP, Exp);
sim = sim/max(abs(sim));

DDP.J = 10000;
[field2, sim2] = pepper(DDP, Exp);
sim2 = sim2/max(abs(sim2));

plot(field, sim)
hold on
plot(field2, sim2)

%% triplet doublet pair
clear all

TDP.J = -6500; % Change the exchange interaction, a value was selected here so that anticrossings appear
TDP.S =[1, 0.5];
TDP.g = [2.003; 2.0059];
TDP.D = [700 -100; 0 0];

TDP.initState = {[0.5 0.2 0.2 0.5 0.1 0.1], 'coupled'}; % switch on and off spin polarisation
TDP.pop_d = [0.2 0.2];
TDP.lw = 0.4;

Exp.Range = [320, 380];
Exp.mwFreq = 9.75;
Exp.Harmonic = 0;
Exp.Temperature = 80;

[field, sim] = pepper(TDP, Exp);
sim = sim/max(abs(sim));

TDP.J = 10000;
[field2, sim2] = pepper(TDP, Exp);
sim2 = sim2/max(abs(sim2));

plot(field, sim)
hold on
plot(field2, sim2)

But I have also observed this behaviour with other parameter sets. Shown here are two extreme examples where the anticrossing is in the observed range and where the spectrum looks clearly strange. However, I have observed that the closer you get to the area with anticrossings, the more ‘wrong’ the spectra become. I used my own programme as a comparison. This calculates the spectra completely differently and more inefficiently than EasySpin. In areas that are far away from anticrossings, the spectra match. The closer you get to the anticrossings, the more differences there are, which are barely noticeable at first and then become increasingly clear.

I hope I have been able to describe the problem clearly enough. If not, please ask me what I mean. I look forward to an answer (even if the error is on my side when setting up the scripts) - thank you in advance!

[Stefan Stoll]

This is a known and unfortunate downside of the (otherwise very accurate and efficient) EasySpin way of simulating powder EPR spectra.

It occurs whenever state characters change rapidly as a function of orientation, for example around anticrossings such as in your cases.

There are actually two separate issues: (1) part of one transition is missed by pepper's transition pre-selection, (2) the orientational interpolation algorithm is inaccurate around anticrossings.

In these cases, one has to tell pepper to (1) include all possible transitions in the simulation and (2) use a lot of orientations. Do this by using

Code: Select all

Opt.Threshold = 0;  % include all transitions (a small number like 1e-9 works as well)
Opt.GridSize = [200 0];  % use a lot of orientations (the larger the better, but the slower as well)
pepper(DDP,Exp,Opt);

[theresia_q]

Thank you very much for the reply and I apologize for reporting a known issue here.

Nevertheless, I have another question: I have implemented your suggested solution and in the case of the boltzmann-populated doublet-doublet pair it works great. However, if I select a singlet precursor, suddenly only one peak comes out:

ddp_2.png (31.48 KiB) Viewed 13314 times

ddp_2_singlet_precursor.png (34.64 KiB) Viewed 13314 times

I can well imagine that I made a mistake here myself, so here again is the code I used to create it. If that's the case, I'm sorry - otherwise I hope it's useful for you to learn about the problem.

Code: Select all

% Spin system
DDP.J = 9500; % Change the exchange interaction, a value was selected here so that anticrossings appear

DDP.S = [1/2 1/2];
DDP.g = [2.0020; 1.9999];
DDP.dip = 350;

DDP.lwpp = 2; 
DDP.initState = 'singlet'; % switch on and off spin polarisation

% Experimental parameters
Exp.Range = [310.777, 374.777];
Exp.mwFreq = 9.6;
Exp.Harmonic = 0;
Exp.Temperature = 80;

Opt.Threshold = 0;
Opt.GridSize = [300, 0];


[field, sim] = pepper(DDP, Exp, Opt);
sim = sim/max(abs(sim));

DDP.J = 10000;
[field2, sim2] = pepper(DDP, Exp);
sim2 = sim2/max(abs(sim2));

plot(field, sim)
hold on
plot(field2, sim2)
hold off

Furthermore, your solution to the problem in the case of the triplet-doublet pair works neither in the Boltzmann-populated case nor in the spin-polarized case:

tdp_2.png (45.93 KiB) Viewed 13314 times

Here I increased the GridSize until the spectrum was no longer changed. The wavy parts from my previous post are now gone. Nevertheless, the image does not look right, e.g. the signal is too wide and it is also not symmetrical. I don't understand where these signals come from and apparently they are not dependent on the number of grid points.

[Matt Krzyaniak]

Hi, perhaps I can help to answer your issues, and no these aren't exactly bugs, but more of a consequence of the spin system and the built in ES assumptions.

So in your first case, you're comparing two simulations with slightly different but large exchange couplings. With a singlet initialization state the singlet level is purely populated, but due to mixing driven by your spin Hamiltonian, you get a small degree of S-T0 mixing.

First, in the second simulation, when J = 10 GHz, if you were to not normalize the spectra, compare the amplitudes of singlet initialization and not (also turn off the temperature), you go from a relative amplitude of 120 to 10^-5 when you turn on the initialization. Although the mixing allows for some T0 character, its not really allowed and is not likely very strong nor polarized.

Second, when you drop the exchange coupling some, in addition to that S-T0 mixing, you get an anti-crossing point and S-T-1 mixing. The sharp spike that you observer comes from one or possibly a range of orientations where there is a significant amount of mixing, giving you a giant absorptive peak. The relative amplitude of that absorptive peak is going to depend on how close to the exact orientation of maximum mixing that is represented in your orientation grid. If you have a super fine grid, >300 pts you'll get a large amplitude, a coarse grid and the amplitude will be relatively small and if its large enough or you know that exact orientation, you can likely get an amplitude comparable to the non-initialized spectrum.

Part of the problems around anticrossing points is that they are treated very discretely in ES, and so it may be difficult to capture exactly what happens in a real system. Inhomogenous hyperfine interactions can likely broaden or even drive the mixing at the anticrossing point, something that is not easy to include in ES. One workaround would be to use a small distribution of J, not a strain or line broadening, but a weighted sum over a range of J.

I can't really help to explain what is happening in your triplet-doublet pair spectra without an example script. The advice I can give, if you don't have strong experimental evidence for an anticrossing( a strong peak where you don't expect one) ignore them and either set J really large(or small) or shift your coupled system to an "infinite J" and turn it into an S=3/2. Also it may be worth playing around with the other polarization modes, rather than just singlet, such as T0 or assigning population to the levels in the Sys.initState = {[x x x x x],'eigen'};

[theresia_q]

Thanks for your detailed answer, Matt! It definitely makes sense to me for the doublet-doublet pair.

Regarding the triplet-doublet pair: The code I used to generate the spectrum is the same as in my first post extended by

Code: Select all

Opt.Threshold = 0;
Opt.GridSize = [300, 0];

I totally understand that it is certainly not an easy problem to solve to generate the spectrum correctly in the case of anticrossings. I have also tested other simulation packages that calculate in a similar way to Easyspin and the problem also occurs there. So it really doesn't seem to be simple...

Nevertheless, I would like to point out again that after numerous further simulations and comparisons, I have now realised that the spectra near the anticrossings look very strange (as in my previous post in the figure) and also show deviating peak positions a good distance away from the anticrossings. Wouldn't that mean that the resonant fields are not calculated correctly in some cases? It is true that the simulations show no problems when using very large Js. However, I am trying to analyse the behaviour of the spectra near the anticrossings (so using large J is not an option for me) and I have noticed that the spectra do not look quite correct even a good distance away from them.

Even if, like me, you don't have the admittedly rare goal of investigating the fundamental behaviour of spectra near the anticrossings, this becomes a problem when investigating systems in the intermediate coupling regime. You don't necessarily have to have an anticrossing there, but J is not yet so large that the eigenstates are well separated. Consequently, you get a spectrum that has the shape of a quartet and seems to look right. However, as described above, this is not entirely the case and this could happen without you realising it.

With regard to the different initial polarisation modes: This is independent of that. The problem also occurs for the ‘ eigen’ precursor.