This user guide explains how to use spidyan for advanced pulse EPR simulations and to investigate spin dynamics.

It is easier to follow this guide, if you are familiar with part I and part II of the guide on saffron.

For the pulse definitions the documentation page of pulse and the corresponding guide are a great starting point.

This guide covers the following topics:

• Running the simulation
• Inversion pulses
• 3D plot of adiabatic passage
• Detection
• Enhancement of the central transition of a high-spin system
• Selective excitation and detection operators
• Rabi oscillation
• Spin echo of a nitroxide
• Spin echo of a Gaussian line
• Simultaneous irradiation at two frequencies
• Tips and Tricks

A spidyan script has the same layout as a saffron script. First a spin system is defined:

Sys.ZeemanFreq = 9.5; % in GHz

After the spin system follows the experiment. Here we are creating a two-pulse sequence:

p90.Flip = pi/2; % flip angle in rad
p90.tp = 0.02; % pulse length in μs

p180.Flip = pi; % flip angle in rad
p180.tp = 0.04; % pulse length in μs

tau = 0.5; % delay in μs
Exp.Sequence = {p90 tau p180 tau}; % pulse sequence
Exp.mwFreq = 9.5; % carrierfrequency

Here we create two monochromatic rectangular pulses with flip angles π/2 and π that are separated by a free evolution period.

Running the above code with

spidyan(Sys,Exp)

will plot the evolution of the electron coherence during the entire sequence (this includes during the pulses).

You might be expecting an echo to form at some point, after all this is the common two-pulse echo sequence. But keep in mind, an echo is a property of a distribution of spins with different resonance frequencies - the echo forms at a certain point of time, when all the spins are refocused. Since we called spidyan with a single resonance frequency, no echo is to be expected.

This section is going to show three different inversion pulses and how to plot their effects on the spin. First, we start again with the spin system.

Sys.ZeemanFreq = 9.5; % in GHz

Exp.DetOperator = 'z1';
Exp.mwFreq = 9.5; % GHz

Since we are interested in the inversion of the spin we chose Ŝ_z as detection operator.

First, let us start with a commong monochromatic rectangular pulse:

Monochromatic.tp = 0.06; % pulse length in μs
Monochromatic.Flip = pi; % flip angle in rad

This is the minimal required input for monochromatic pulse. If no Pulse.Type is given, a rectangular pulse is assumed by default. The same is true for the parameter Pulse.Frequency, which is set to 0 by default, corresponding to a pulse with the frequency of the carrier frequency Exp.mwFreq.

For completeness, here is an equivalent pulse definition:

Monochromatic.Type = 'rectangular';
Monochromatic.tp = 0.06;
Monochromatic.Flip = pi;
Monochromatic.Frequency = 0;

Next, let us build Exp.Sequence, run the simulation and then plot it

Exp.Sequence = {Monochromatic};

[TimeAxis, Signal] = spidyan(Sys,Exp);

% plotting: 
figure(1)
clf
hold on
xlabel('t (ns)')
ylabel('< S_z>')
axis tight
ylim([-1 1])

plot(TimeAxis*1000,real(Signal));

Next, let us look at what happens to a spin during adiabatic passage. We define the pulse and overwrite current Exp.Sequence

LinearChirp.Type = 'quartersin/linear'; 
LinearChirp.trise = 0.030;  % rise time μs
LinearChirp.tp = 0.2; % pulse length μs
LinearChirp.Frequency = [-50 50]; % frequency band in MHz
LinearChirp.Flip = pi; % flip angle in rad

% Update sequence
Exp.Sequence = {LinearChirp};

According to pulse, this gives a chirp pulse with a length of 200 ns and a bandwidth of of 100 MHz, centered around the carrier frequency Exp.mwFreq. To ensure a well-defined excitation profile in the frequency domain, the edges in the time domain are smoothed with a quarter sine, which has a length of 30 ns.

Run and plot this with

[TimeAxis, Signal] = spidyan(Sys,Exp);

% Plotting
plot(TimeAxis*1000,real(Signal));

and you will see how the spin follows the effective field during adiabatic inversion.

Once again, we define a new pulse and plot the simulation:

HS.Type = 'sech/tanh';
HS.beta = 10; % truncation parameter of 'sech'
HS.n = 1; % exponent of the secant function argument
HS.tp = 0.2; % pulse length
HS.Frequency = [-50 50]; % excitation band, GHz
HS.Flip = pi; % flip angle in rad

% Update sequence
Exp.Sequence = {HS};

[TimeAxis, Signal] = spidyan(Sys,Exp);

% Plotting
plot(TimeAxis*1000,real(Signal));

Although trajectory is slightly different, both adiabatic pulses end up in the same spin state.

The quality, of adiabatic pulses is in general described by a parameter called critical adiabaticity Qcrit, which measures the ability of the pulse to invert the spin and depends on the amplitude and the sweep rate. The flip angle and Qcrit are related through the Landau-Zener-Stückelberg-Majorana equation:, e.g, a Qcrit of 5 corresponds to a π pulse and a Qcrit of about 0.7 to a π/2 pulse. Hence, it is not only possible to define the pulses with Pulse.Flip, but instead you can provide Pulse.Amplitude or Pulse.Qcrit for frequency-swept pulses.

Adiabatic pulses have the big advantage, that they invert all spins within their excitation band equally. You can see this, by changing the resonance frequency of the spin Sys.ZeemanFreq from 9.500 by to 25 MHz to Sys.ZeemanFreq = 9.525 and then comparing the three pulses.

A second advantage of adiabatic pulses is, that overturning of the spin is not possible, something that is common with rectangular pulses, see this example. This means, no matter how strong the amplitude of your pulse gets, the flip angle can not became larger than π. You can try this by removing Par.Flip from the pulse definitions, and instead using Par.Amplitude.

The aim of this example is to create a three-dimensional depiction of the spin trajectory during adiabatic passage. We start with the spin system:

Sys.ZeemanFreq = 9.500; % resonance frequency of the spin in GHz

Next, the pulse and experiment definition:

Pulse.Type = 'quartersin/linear';
Pulse.trise = 0.015; % rise time in μs
Pulse.Qcrit = 5;  % critical adiabaticity
Pulse.tp = 0.2; % pulse length in μs
Pulse.Frequency = [-100 100]; % frequency band in MHz

Exp.Sequence = {Pulse}; % μs
Exp.mwFreq = 9.500; % GHz

For the three-dimensional plot we need the expectation values of Ŝ_x, Ŝ_y and Ŝ_z. Since Ŝ+ = Ŝx + i Ŝy we can use this to our advantag and use for the detection:

Exp.DetOperator = {'z1','+1'}; % detection operators
Exp.DetFreq = [0 9.5]; % downconversion frequency in GHz

Exp.DetFreq contains the frequencies for downconversion in GHz. Since Ŝ_z does not contain any rotating components, the down conversion has to be switched off for it, by setting it to 0.

We can run the simulation and plot it: [TimeAxis, Signal] = spidyan(Sys,Exp);

% Plotting
figure(1)
clf
plot(TimeAxis*1000,real(Signal));
xlabel('t (ns)')
axis tight
ylim([-1 1])
ylabel('< S_i>')
legend(Exp.DetOperator)

In order to get a three-dimensional plot, we need use the real and imaginary parts of <Ŝy>:

figure(2)
clf
plot3(real(Signal(:,2)),imag(Signal(:,2)),real(Signal(:,1)));
xlabel('< S_x >')
ylabel('< S_y>')
zlabel('< S_z>')
grid on

More on the theory of adiabatic passage can be found in

Gunnar Jeschke, Stephan Pribitzer, Andrin Doll
Coherence Transfer by Passage Pulses in Electron Paramagnetic Resonance Spectroscopy
The Journal of Physical Chemistry B 2015, DOI: 10.1021/acs.jpcb.5b02964

and

Philipp Spindler, et al.
Shaped Pulses in EPR
eMagRes 2016, DOI: 10.1002/9780470034590.emrstm1520

The use of the detection window Exp.DetWindow was already discussed in part II of the guide on saffron.

With spidyan you get more control over what is detected and when.

The first topic concern detection operators. By default, if Exp.DetOperator is omitted, Ŝ₊ of all electrons is detected.

But using Exp.DetOperator it is not only possible to use several operators at once, you can even specify spin and transition selective operators, using the extended sop syntax. The detection operators are specified in a cell array, e. g.

Exp.DetOperator = {'z1' 'x1' '+2'};

uses Ŝ_z and Ŝ_x of the first spin of the spin system and Ŝ₊ of the second spin.

For each of the components of Exp.DetOperator a frequency for downconversion in GHz must be provided. Detection operators that correspond only to elements on the diagonal of the detection operator do not need to be down converted. Such operators are any Ŝ_z-operators, written like 'z1', polarization operators between two levels, specified for example with z(1|2)1 and population operators selective operators, specified with like e(2).

Assuming that both spins are of the same type, we could write for the above example:

Exp.DetFreq = [0 9.5 9.5]; % down conversion frequency in GHz

This corresponds to no down conversion of 'z1' and down conversion with 9.5 GHz for 'x1' and 9.5 GHz for 'x2'.

For counter-rotating operators, e.g. Ŝ_-, the sign of the corresponding element in DetFreq has to reflect this:

Exp.DetOperator = {'+1' '-1'}
Exp.DetFreq = [9.5 -9.5]; % down conversion frequency in GHz

Instead of using Exp.DetWindow, which only allows detection after the last pulse, you can also use Exp.DetSequence, which allows you to detect during selected events.

If neither Exp.DetWindow nor Exp.DetSequence are defined, spidyan detects the entire sequence by default.

Detection is controlled through booleans for each element in Exp.Sequence separately. Here are a few examples:

Exp.Sequence = {p90 tau p180 2*tau}; % the pulse experiment

Exp.DetSequence = 1 % detect the entire sequence
Exp.DetSequence = 0 % no detection at all
Exp.DetSequence = [0 0 0 1] % only detect during the last event in Exp.Sequence
Exp.DetSequence = [0 1 0 1] % detection during the second and fourth element in Exp.Sequence
Exp.DetSequence = [0 0 1 0] % detection during the third event

Keep in mind, that detection not only requires a significant amount of computing time, it also keeps spidyan from using some of it speed-up tricks. Hence, if you run big simulation with a lot indirect dimensions and data points, you might want to think about reducing detection to a minimum for the sake of speed.

This example is going over the enhancement of the central transition of a S = 3/2 system for a single spin packet. For this we are going to use two frequency swept pulses that will transfer magnetization from the outer transitions to the central transition. We can follow the polarization enhancement by using transition selective operators.

First, let us define a spin system with some zero-field splitting:

Sys.S = 3/2;
Sys.ZeemanFreq = 33.500; % GHz
Sys.D = 166;  % MHz

Next we set up two chirp pulses that are being swept from the outside of the spectrum towards the center and put them right after each other in Exp.Sequence:

Pulse1.Type = 'quartersin/linear'; % make it a chirp
Pulse1.trise = 0.05; % smooth edges in time domain with a quarter sine, in mus
Pulse1.Qcrit = 10; % use critical adiabaticity instead of Par.Flip or Par.Amplitude
Pulse1.tp = 1; % pulse length in mus
Pulse1.Frequency = [500 170]; % frequency band, relative to Exp.mwFreq

Pulse2.Type = 'quartersin/linear'; % make it a chirp
Pulse2.trise = 0.05; % smooth edges in time domain with a quarter sine, in mus
Pulse2.Qcrit = 10; % use critical adiabaticity instead of Par.Flip or Par.Amplitude
Pulse2.tp = 1; % pulse length in mus
Pulse2.Frequency = [-500 -170]; % frequency band, relative to Exp.mwFreq

% Sequence
Exp.Sequence = {Pulse1 Pulse2};

We then set the carrier frequency and specify our detection operators:

Exp.mwFreq = 33.5; % GHz

% The detection operators detect polarization between (1) levels 1 and 2,
% levels 2 and (3) levels 3 and 4
Exp.DetOperator = {'z(1|2)' 'z(2|3)' 'z(3|4)'};

We make use of the transition-selective operators from sop: With 'z(1|2)' polarization between |3/2> and |1/2> is detected, with 'z(2|3)' the central transition (|1/2> and |-1/2>) and with 'z(3|4)' the polarization of the states |-1/2> and |-3/2>.

Running and plotting it with

[TimeAxis, Signal] = spidyan(Sys,Exp);

% plotting
figure(1)
clf
plot(TimeAxis,-real(Signal));
xlabel('t ({\mu}s)')
axis tight
ylim([-1 3])
ylabel('< S_i>')
legend(Exp.DetOperator)

shows how the first pulse shuffles polarization from one of the outer transitions to the central transition, and increasing the polarization on the central transition from 1 to to 2. The second pulse does the same with the other out transition, raising the polarization of the central transition to 3! Try this with a S = 5/2 or S = 7/2 and see if you can achieve even larger enhancements - you might have to adapt the excitation bands of your pulses and add the correct detection operators!

Now, this simulation was done for a single spin packet (e.g. a single crystal) only. In reality, you encounter a distribution of spins (the spectrum) paired with a distribution of zero-field splitting values which make it harder to enhance the central transition. A real life example with Gd(III) can be found in

Andrin Doll, et al.
Sensitivity enhancement by population transfer in Gd(III) spin labels
Phys. Chem. Chem. Phys. 2015, DOI: 10.1039/C4CP05893C

This example explores how spin selective excitation and detection operators are used. First we start with the specification of a spin system. Since we are interested in seeing the pulses and not spin dynamics, we are going to set the spin-spin coupling to zero.

Sys.S = [1/2 1/2];
Sys.ZeemanFreq = [9.500 9.500]; % GHz

Sys.J = 0; % spin-spin coupling, MHz

Next we are going to define pulses and then put them into a pulse sequence:

P90.tp = 0.016; % mus
P90.Flip = pi/2; % rad

P180.tp = 0.032; % mus
P180.Flip = pi; % rad

tau = 0.5; % mus

% Sequence
Exp.Sequence = {P90 tau P180 tau P180 tau};

Exp.mwFreq = 9.5; % GHz

To see how the pulses flip our spins, we are going to use selective Ŝ_z operators for both spins:

Exp.DetOperator = {'z1' 'z2'};

Next, let us use some selective excitation operators with Opt.ExcOperator. The indexing in Opt.ExcOperator. corresponds to the pulses, e.g. the second element in Opt.ExcOperator concerns only the second pulse. For declaring special excitation operators, we can use the usual sop syntax.

Opt.ExcOperator = {'x1' [] 'x2'}; 

Now, the first pulse will use Ŝ_x for the first spin, the second pulse excites all spins in the system and the last pulse only excites the second spin. Run and plot it with

[TimeAxis, Signal] = spidyan(Sys,Exp,Opt);

% plotting
figure(1)
clf
plot(TimeAxis,real(Signal));
xlabel('t ({\mu}s)')
axis tight
ylim([-1 1])
ylabel('< S_i>')
legend(Exp.DetOperator)

and you will see the spin flips.

Of course you can use transition selective excitation operators by using the sop syntax. If it is not possible to build the excitation operator with sop, it is even possible to manually design the excitation operator in matrix form, and then use that as input:

Sys.S = 3/2;

% transition selective operator using sop syntax:
Opt.ExcOperator = {'x(2|3)'}; 

 % manually define x(2|3) for S = 3/2:
Op = [0 0   0   0;   
      0 0   0.5 0;
      0 0.5 0   0;
      0 0   0   0];

Opt.ExcOperator = {Op}; 

This is example is going how to setup a rabi oscillation with spidyan using indirect dimensions. For more info on indirect dimension, please refer to the documentation. Once again, we start with the spin system: and then set up a very short 1 ns rectangular pulse.

% Spin System
Sys.ZeemanFreq = 9.500; % GHz

In order to observe Rabi oscillations, we want a rectangular pulse with a constant amplitude. We then are going to use the indirect dimension to increase the length of the pulse and detect the outcome of the experiment. Let us start with the pulse and Exp.Sequence

% experiment setup
Pulse.Type = 'rectangular';
Pulse.tp = 0.001; % pulse length in μs
Pulse.Amplitude = 30; % pulse amplitude in MHz

Exp.Sequence = {Pulse}; 
Exp.mwFreq = 9.5; % GHz

Next, detection. As detection operator we use

Exp.DetOperator = {'z1'};

and we are interested in the expectation value right after the pulse. The easiest way is to use Exp.DetWindow:

Exp.DetWindow = 0; % position of the detection window in μs

Since we are not providing an interval for the detection window, we are going to do single point detection. The value 0 corresponds to the time in μs after the end of the last element in Exp.Sequence. For more information on how to use Exp.DetWindow, see the documentation.

Alternatively, we could have adapted Exp.Sequence and used Exp.DetSequence (see here and here):

Exp.Sequence = {Pulse 0}; % append an event with length 0 μs

Exp.DetSequence = [0 1]; % detect only during the the second element in Exp.Sequence

The indirect dimension, where we change the pulse length, can be added through:

Exp.nPoints = 99 % number of data points 
Exp.Dim1 = {'p1.tp' 0.001}; % increment length of the pulse by 1 ns per step

With this, we are going to record 99 data points, with the pulse length ranging from 1 ns to 100 ns:

spidyan(Sys,Exp);
ylim([-1 1]) % to ensure right limits in the plot

Now you can investigate how the frequency of the Rabi oscillation depends on the frequency difference between the microwave pulse and the resonance frequency of the spin by changing Sys.ZeemanFreq.

Or you can add some relaxation by adding:

Sys.T1 = 2; % longitudinal relaxation time in μs
Sys.T2 = 1; % transverse relaxation time in μs

Opt.Relaxation = 1;

% then run
spidyan(Sys,Exp,Opt);

Per call spidyan can process a single spin packet. If you want to simulate a pulse sequence with an entire spectrum, you would usually use saffron, unless you want to use special excitation operators, change the detection operator or detect between and during pulses. In such a case you can manually define an orientation loop which calls spidyan.

We start again with the spin system, but this time we are using Sys.g instead of Sys.ZeemanFreq. Since we are going to manually rotate the interactions, we are going to set them up as tensors:

% Spin System
Sys.S = 1/2;
Sys.g = diag([2.00906 2.0062 2.0023]);
Sys.A = diag([11.5 11.5 95]);
Sys.Nucs = '14N';

We then determine the grid with the EasySpin functions hamsymm and sphgrid.

Symmetry = hamsymm(Sys);
[phi,theta,Weights] = sphgrid(Symmetry,20);

Then we set up frequency-swept pulses, that cover the entire spectrum:

Chirp90.Type = 'quartersin/linear'; 
Chirp90.trise = 0.030;           % edge smoothing, μs
Chirp90.tp = 0.200;              % pulse length, μs
Chirp90.Flip = pi/2;             % flip angle, rad
Chirp90.Frequency = [-120 120];  % excitation band, GHz

Chirp180.Type = 'quartersin/linear';
Chirp180.trise = 0.030;             % smoothed edges, μs
Chirp180.tp = 0.100;                % pulse length, μs
Chirp180.Flip = pi;                 % flip angle, rad
Chirp180.Frequency = [-120 120];    % excitation band, GHz

tau = 0.5;  % μs
Exp.Sequence = {Chirp90 tau Chirp180 tau+Chirp180.tp}; 

Exp.mwFreq = 9.1;   % GHz 
Exp.Field = 324.9;  % mT

% Detection:
Exp.DetWindow = [-0.05 0.05];  % μs
Exp.DetOperator = {'+1'};
Exp.DetFreq = 9.1;

If you are not sure about where to place your pulses, keep in mind that you can call pepper with a Sys and Exp.Field, which is going to plot the spectrum for you. For detection we are detecting the electron coherence in a 0.1 μs window centered around the end of the last element in Exp.Sequence. Since we are using frequency-swept pulses, the position of the echo is not at τ after the end of the π pulse, but we also need to take the pulse length into account.

With this we can loop over the orientations and accumulate the signal:

Signal = 0; % initialize Signal
for iOrientation = 1 : numel(Weights)
  
  Sys_ = Sys; % create temporary copy of Sys
  
  R = erot(phi(iOrientation),theta(iOrientation),0); % rotation matrix
  
  Sys_.g = R'*Sys_.g*R; % rotate Sys.g
  Sys_.A = R'*Sys_.A*R; % rotate Sys.A
  
  [t, signal_] = spidyan(Sys_,Exp);
  
  Signal = Signal + signal_*Weights(iOrientation); % sum up signals
end

After completion, the signal can be plotted:

% plotting
figure(1)
clf
plot(t,real(Signal))
xlabel('Transient (\mus)')
ylabel('Signal (arb.u.)')

If you are investigating spin dynamics, it is often not necessary to simulate entire spectra. Instead you can use any distribution of spin packets (a distribution of different resonance frequencies is needed for an echo to form).

Here we are going to cover how to loop spidyan over a Gaussian line. The resonance frequencies of the spins are centered around 9.5 GHz, with a standard deviation of 10 MHz. Spins are sampled from 9.450 to 9.550 GHz, with a sampling step 0.5 MHz, resulting in a total of 201 spin packets.

CenterFrequency = 9.5; % center frequency of Gaussian distribution, GHz
GWidth = 0.01;     % width of Gaussian distribution, GHz
FreqStart = 9.45;  % starting value for sampling
FreqEnd = 9.55;  % final value for sampling
Sampling = 0.0005;   % stepsize for sampling
ZeemanFreqVec = FreqStart:Sampling:FreqEnd; % vector with resonance frequencies
P = exp(-((CenterFrequency-ZeemanFreqVec)/GWidth).^2); % probabilities
P = P/trapz(P); % normalization
nSpinpackets = length(ZeemanFreqVec);

For each of the spin packets we now have a resonance frequency and a probability.

After setting up the pulse sequence

Pulse90.Type = 'rectangular';
Pulse90.tp = 0.025; % pulse length, mus
Pulse90.Flip = pi/2; % flip angle, rad

Pulse180.Type = 'rectangular';
Pulse180.tp = 0.025; % pulse length, mus
Pulse180.Flip = pi/2;  % flip angle, rad

Exp.Sequence = {Pulse90 0.25 Pulse90 0.5}; 
Exp.mwFreq = 9.5; % GHz

Exp.DetSequence = [0 0 0 1]; % detect the last event

Exp.DetOperator = {'+1'}; % detect electron coherence
Exp.DetFreq = 9.5; % GHz

Now, all that is left, is looping over the spins. In contrast to the example with the nitroxide, we do not need to rotate any tensors anymore, and can just load the resonance frequency and then plot the accumulated signal:

Signal = 0; % initialize the signal

% Loop over the spin packets and sum up the traces ------------
for i = 1 : nSpinpackets
  
  Sys_.ZeemanFreq = ZeemanFreqVec(i); % Set Zeeman frequency
  
  [t, signal] = spidyan(Sys_,Exp);
  
  Signal = Signal + signal*P(i);
  
end

figure(1)
clf
plot(t,abs(Signal));
xlabel('t (\mus)')
axis tight
ylim([0 1])detect electron coherence
Exp.DetFreq = 9.5; % GHz

Simultaneous irradiation (e.g. cross-polarization or decoupling) at two different frequencies is not implemented at this point, but can still be achieved: The trick is to use Exp.mwPolarization in combination with custom excitation operators. This allows us to assign the real part of the wave to an excitation operator that acts on one spin or subset of spins (frequency) and the imaginary part to a second spin or subset of spins (frequency). Look at the following example:

  CombinedPulses.IQ = real(mwPulse) + imag(rfPulse);
  CombinedPulses.t = t_mwPulse);

  Opt.ExcOperator{1} = sop([1/2 1/2],'xe') + 1i*sop([1/2 1/2],'ex');
  
  Exp.mwPolarization = 'circular'; 
  
  Exp.Sequence = {CombinedPulses};
  
  spidyan(Sys,Exp,Opt)

In the next step we assign spin selective operators: The electron spin is to get excited only real part of the pulse and the nucleus by the imaginary part. With Exp.mwPolarization = 'circular' we are letting spidyan (or saffron) know that it should expect a wave form that contains a real and an imaginary part.

If the above code is now run, the simulation will use the real part of the wave form to propagate the electron spin and the imaginary part to propagate the nucleus, as if we were irradiating the spin system with two pulses at different frequencies at the same time.

If you are interested in the expectation values of Ŝ_x, it is usually beneficial to use Ŝ₊ as the actual detection operator and then take the real part of the obtained signal. This removes artifacts at the beginning and end of the time traces that are introduced when translating a purely real signal during the signal processing.

Ŝ₊ and Ŝ_- are not Hermitian operators. Hence, if you use Ŝ₊ as detection operator, the signal that is returned will in fact correspond to ⟨Ŝ_-⟩. The real part of the signal (Ŝ_x) is nof affected by this. But if your data processing involves the imaginary part of the signal and you encounter frequencies with a wrong sign, you might want to consider using Ŝ_- as detection operator instead.

In some simulations you may observe echoes in your timetrace that should not be there. Besides physical reasons (incomplete phase cycling, different refocusing conditions,...), such artifact echoes can arise from an aliasing effect if the spectrum (e.g. the nitroxide spectrum in the example above) is insufficiently sampled. This can easily be checked: By increasing the number of samples from your distribution the artifact echoes should move, spread further apart or vanish, while physical echoes stay in place. You can also compute the required increment for the spin packet distribution in frequency domain by an inverse Nyquist criterion. If you want to observe up to a time t_max after the first excitation event (first pulse), the frequency increment should not be larger than 1/(2 t_max). Otherwise use Monte Carlo sampling of spin packet frequencies, which will convert the artifacts to 'Monte Carlo noise' that averages with increasing number of trials.

If the signal processing fails, spidyan returns the signal in the simulation frame. This can happen if you try to translate Ŝ_z or provide a wrong frequency. Instead of adapting down conversion and rerunning the entire simulation, you can call signalprocessing with the returned signal and the correct down conversion frequencies. More info here.