spidyan(Sys,Exp)
spidyan(Sys,Exp,Opt)
sig = spidyan(...)
[t,sig] = spidyan(...)
[t,sig,out] = spidyan(...)

Array with Zeeman frequencies (in GHz) of electron spins defined in S of the spin system structure. It is possible to use ZeemanFreq in combination with g. If a non-zero Zeeman frequency is found, any g-values given for this electron spin are ignored.

Longitudinal relaxation time in microseconds. If it is a scalar, the relaxation time is applied to all transitions. Transition specific relaxation times can be provided in the form of a matrix. See section Relaxation times below for more details. Relaxation times for transitions that are not defined (or 0) are automatically set to 10¹⁰ microseconds. In order for spidyan to simulate relaxation effects, it is necessary to switch relaxation on with Opt.Relaxation (see below).

Transverse relaxation time in microseconds. If it is a scalar, the relaxation time is applied to all transitions. Transition specific relaxation times can be provided in the form of a matrix. See section Relaxation times more details. Relaxation times for transitions that are not defined (or 0) are automatically set to 10¹⁰ microseconds. In order for spidyan to simulate relaxation effects, it is necessary to switch relaxation on with Opt.Relaxation (see below).

It is possible to run spidyan simulations with a custom initial state, which must be provided in matrix form. By default, the initial state assumes the high temperature approximation for all electron spins. The structure of density matrices (ordering of operators in terms of product Zeeman basis states) is explained at sop.

Sys.S = 1/2;				% Creates an isolated spin 1/2
Sys.initState = [0 1/2; 1/2 0]; 	% Start simulation from Sx

For simulations with relaxation it is possible to provide an equilibrium state eqState that the system relaxes to. This state has to be a density operator matrix form. By default the initial state is used as equilibrium state. The equilibrium state is only used if relaxation is active, else it is ignored.

Magnetic field (in mT) at which the experiment is performed, needs to be provided only if the spin system is defined with Sys.g. If Sys.ZeemanFreq is used, Field can be omitted.

A cell array that contains the pulse definition and the inter-pulse delays. Available pulse types and their parameters are described at pulse. For a pulse it is necessary to insert a structure containing the pulse definition, free evolution events are declared with a scalar that corresponds to their length in microsecond:

P90.Type = 'rectangular';  	% Define Pulse, see examples for more details
P90.tp = 0.032;		  	% Pulse Length in microseconds
P90.Flip = pi/2;		  	% Flip angle of the pulse

P180.Type = 'rectangular';  	% Define Pulse, see examples for more details
P180.tp = 0.032;		% Pulse Length in microseconds
P180.Flip = pi;		  	% Flip angle of the pulse

Exp.Sequence = {P90 0.2 P180 0.4}; % Creates a two pulse (echo) sequence, with an inter-pulse 
				   % delay of 0.2 microseconds and a final free-evolution period 
				   % of 0.4 microseconds

The delays are defined to go from the end of one previous element in Sequence to the beginning of the next (unlike in Bruker spectrometers).

EPR spectrometer frequency in GHz. All frequencies in the pulse definition (e.g. Pulse.Frequency) need to be defined relative to it:

Exp.mwFreq = 33.5;

Pulse.Frequency = [-250 250];  % 500 MHz sweep

Cell array that contains the individual phase cycles, with the indexing corresponding to the pulse indices. The phase cycle itself needs to be given in array form, where rows correspond to the phase cycling steps. In each row the first element is the phase, and the second element the detection phase.

PC = [0, 1; pi, -1];   		% phase cycle
Exp.PhaseCycle = {[] PC};	% phase cycles the second pulse

Phase cycling that steps more than one pulse simultaneously is not available through the PhaseCycle structure, but can be achieved through adding an additional indirect dimension to the experiment (see Dim).

Optional, can be 'linear' or 'circular', default is 'linear'. If mwPolarization is 'circular', the excitation operator also uses the imaginary part of the waveform. The default excitation then becomes real(IQ) Ŝ_x + imag(IQ) Ŝ_y.

Exp.mwPolarization = 'circular';     % mwPolarziation set to circular for all pulses

mwPolarization can also be used in combination with custom excitation operators (Opt.ExcOperator), if they have a imaginary component. Keep in mind, that using 'circular' does not allow spidyan to use some of its speed up tricks and can lead to significantly longer simulation times.

The field DetWindow controls detection and is given in microseconds: The zero-time of the detection is taken as the end of the last element in Exp.Sequence.

To do single point detection, it is sufficient to tell DetWindow the time of the detection point. E.g. to detect at the center of the echo after a two pulse experiment:

tau = 0.5
Exp.Sequence = {pulse90 tau pulse180 tau};   	% a two-pulse echo

Exp.DetWindow = 0;	% detect a single point at the end of the second delay (center of the echo)

The detection position can be set to anywhere at the end of the sequence, as long as it does not overlap with a pulse:

Exp.DetWindow =  0.2;	% detect a single point 200 ns after the second delay
Exp.DetWindow = -0.2;	% detect a single point 200 ns before the end of the second delay

In order to detect a transient, DetWindow has to be provided with a start and end value:

tau = 0.5
Exp.Sequence = {pulse90 tau pulse180 tau};   	% a two-pulse echo

Exp.DetWindow = [-0.1 0.1];	% detect a from -100 ns to 100 ns centered around of the end of tau

Of course, it is also possible to detect a transient anywhere after the last pulse, as long is it does not overlap with it:

Exp.DetWindow = [0 0.3];	% start detection at the end of the last event and detect for 300 ns
Exp.DetWindow = [-0.3 -0.1];	% start detection 300 ns before the end of the last event and detect for 200 ns

If you want to detect an FID or have a fixed experiment length (e.g. 3p ESEEM) you can also write:

Exp.Sequence = {pulse90 0.5 pulse90 0.8 pulse90};   	% a three-pulse ESEEM sequence (the last element is a pulse)

Exp.DetWindow = 0.5;	% detect a single point after 500ns after the last pi/2 pulse (center of echo)

Exp.DetWindow = [0.250 0.750] % detect a transient, centered around the echo

Exp.DetWindow = [0 0.100]    % detect an FID after the last pulse

A vector for advanced control of the detection, e.g. detect during pulses. Detection can be set for the entire simulation (meaning detection is on for all events, even pulses) or for specific events:

Exp.DetSequence = true;   	% the entire sequence is detected
Exp.DetSequence = [0 0 1 1];	% expectation values are only returned the third and fourth element in Exp.Sequence
Exp.DetSequence = 0;		% no detection

Besides booleans, Exp.DetSequence also accepts character arrays:

Exp.DetSequence = 'all';   	% the entire sequence is detected
Exp.DetSequence = 'last';	% the last element in Exp.Sequence is detected

If the field DetSequence is not provided, detection is switched off by default. Detection operators can be defined in Exp.DetOperator. If DetSequence is active, but no detection operator is defined, Ŝ₊ is used for all electron spins. The length of DetSequence has to be either 1 or the same length as Exp.Sequence.

If Exp.DetWindow is provided, any input from DetSequence is ignored.

The field DetOperator is a cell array that contains detection operators. If no detection operator is defined, but detection is active, Ŝ₊ is used for all electron spins. Detection operators can be defined by using the same syntax as in sop (see there):

Exp.DetOperator = {'+1' 'z1'};   % detects Ŝ+ and Ŝz of the first electron spin

Detection operators that can not be defined using the sop syntax, can be provided in matrix form:

Exp.DetOperator = {[0 1; 0 0] [1/2 0; 0 -1/2]};   % the same operators as a above for S = 1/2 in matrix form

For convenience it is possible write a single detection operator without the curly brackets:

Exp.DetOperator = '+1'; % same as {'+1'}

Exp.DetOperator = [1/2 0; 0 -1/2];   % identical to {[1/2 0; 0 -1/2]}

If you are interested in the expectation values of Ŝ_x, it is usually beneficial to use Ŝ₊ or Ŝ_- as 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 Ŝ_- instead.

The signal is detected in the simulation frame, which causes all signals from detection operators that have non-zero off-diagonal elements to have an oscillating component. By translating/shifting the frequency it is possible to center the signal around 0 (baseband) or any other desired frequency. The frequency shift has to given as absolute frequency (even if used in conjunction with FrameShift) in GHz and for each detection operator separately. Indexing in DetFreq corresponds to the ordering in DetOperator:

Exp.DetOperator = {'+1' 'z1'};   % detects  Ŝ+ and  Ŝz of the first electron spin
Exp.DetFreq = [33.5 0]; % shifts Ŝ+ by -33.5 GHz
				 % Ŝz does not contain a rotating component and does not need to shifted

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

Exp.DetOperator = {'+1' '-1'};	% detects Ŝ+ and Ŝ- of the first electron spin
Exp.DetFreq = [33.5 -33.5]; 	% shifts Ŝ+ by -33.5 GHz and Ŝ- by 33.5 GHz

Phase of detection in radians, optional. Has to be specified for each detection operator separately:

Exp.DetOperator = {'z1' '+1'};
Exp.DetPhase = [0 pi]; 		% change detection phase of second detection operator by π

The default value is 0. The phase provided in DetPhase is added to the signal, regardless the type of operator.

A vector that contains the number of points in each dimension, e. g. [10 150] corresponds to 10 points in the first and 150 points in the second dimension.

Exp.nPoints = [10 150]; 	% 10 points along the first indirect dimension, 150 along the second

Dim1, Dim2,... provide the fields that are to be changed along the indirect dimensions (field nPoints). The first data point always uses the values defined initially in the experiment definition. All fields that appear in the pulse definition can be changed, e.g:

Exp.Dim1 = {'p2.Flip' pi/8};        % increments the flip angle of the second pulse by pi/8 each step

For free evolution events only the length can be changed:

Exp.Dim1 = {'d3' -0.1};             % decrements the length of the third delay by 100 ns each step

Several parameters can be changed in one dimension:

Exp.Dim1 = {'p2.Flip' pi/8; 'd3' -0.1}; 	% changes flip angle of pulse and duration of free evolution
Exp.Dim1 = {'p2.Flip,p3.Flip' pi/8};  		% flip angles of 2nd and 3rd pulse are simultaneously stepped  

For experiments that involve one or several moving pulses, the identifier Position can be used. This is only possible for pulses that are not the first or last event in the sequence. Pulses are allowed to cross, but must not overlap.

Exp.Dim2 = {'p2.Position' 0.1};   % moves the second pulse 100 ns back each step in the 2nd dimension

A list of increments can be used by providing vectors with the precomputed increments. All increments are always applied to the initial value of the field (as defined in the Exp) and are not related to each other. Hence, if the value in the experiment definition is desired as first data point of the indirect dimension, the fist element has to be zero:

Exp.nPoints = 4;

% this
Exp.Dim1 = {'d2' [0 0.1 0.2 0.3]};
% is equal to: 
Exp.Dim1 = {'d2' 0.1}; 

Complete freedom is given when it comes to providing the list of increments.

Exp.nPoints = 5;
Exp.Dim1 = {'d2' [0 0.1 0.3 0.65 -0.2]};   % increment of the second delay

However, the program checks that changing lengths of events do not lead to overlapping pulses at any acquisition point.

A special case is the field Par.Frequency in the pulse definition of a frequency-swept or of a monochromatic pulse that is defined by identical initial and final frequency. If one value is provided in Exp.DimX, this is going to be added to both elements in the Par.Frequency field - the pulse is moved in the frequency domain:

Exp.Dim1 = {'p1.Frequency' 5};   % changes both elements in Par.Frequency by 5 MHz per step

An equivalent input as the above is the following, where each the increment for each field in Par.Frequency is given:

Exp.Dim1 = {'p1.Frequency' [5 5]};   % changes both elements in Par.Frequency by 5 MHz per step

This syntax can also be used to change the sweep width:

Exp.Dim1 = {'p1.Frequency' [-5 5]};   			% increases sweep width by 10 MHz each step
Exp.Dim1 = {'p1.Frequency' [0 5]};    			% increment only the final frequency of the frequency sweep

A list of increments can be provided as well, by using ';' as separator:

Exp.Dim1 = {'p1.Frequency' [0 0; 10 10; 30 30; 80 80]};  % vector increment for shifting the frequency range of the pulse
Exp.Dim1 = {'p1.Frequency' [0 0; -10 0; 0 30; -80 80]};  % vector increment changing the excitation range of the pulse

For other pulse parameters that are defined by a vector (e.g., the order of an asymmetric hyperbolic secant (HS) pulse Pulse.n or the list of relatives amplitudes of a Gaussian cascade Pulse.A0), selected elements can be incremented by adding an index to the field name:

Exp.Dim1 = {'p1.A0(3)' 0.1; 'p1.A0(4)' -0.1};   % changes relative amplitudes of the third and fourth pulse in a Gaussian cascade
Exp.Dim1 = {'p1.n(2)' 2};    % increases order of the falling flank of a hyperbolic secant pulse by 2 each step

% Also possible for 'Frequency'
Exp.Dim1 = {'p1.Frequency(2)' 5};    % changes only the final frequency of the pulse
Exp.Dim1 = {'p1.Frequency' [0 5]};   % identical to the above and not a list of increments!

Adding an indirect dimension also allows for simultaneous phase cycling of two or more pulses (something that can not be achieved through the Exp.PhaseCycle structure):

Exp.nPoints = 4;
Exp.Dim1 = {'p2.Phase,p3.Phase' pi/4};   % changes the phase of the 2nd and 3rd pulse by pi/4 each step

In the above example, the output of spidyan will contain the individual transients from each phase cycling step and manual merging of the dimensions (with proper detection phase/sign) is required to obtain the phase-cycled signal.

This is an example for a two-dimensional experiment (e.g. for a two-pulse echo) to optimize two pulse parameters:

Exp.nPoints = [20 15];
Exp.Dim1 = {'p1.Frequency,p2.Frequency' [-2.5 2.5]}; % increase excitation band by 5 MHz each step of both pulses
Exp.Dim1 = {'p1.tp' 2; 'p2.tp' 1}; % step pulse lengths, in µs

Exp.ResonatorFrequency = 33.5; 	% resonator center frequency in GHz
Exp.ResonatorQL = 300;	  	% loaded Q-value

Exp.FrequencyResponse = [Frequency; TransferFunction] 	% transfer function and its correspond frequency axis

Center frequency of resonator in GHz.

The loaded Q-value.

This gives the frequency response of the resonator, in the form Par.FrequencyResponse = [Frequency; TransferFunction] with the (possibly experimental) resonator transfer function in TransferFunction and the corresponding frequency axis in Frequency (in GHz). A complex transfer function input in Par.FrequencyResponse is used directly in the bandwidth compensation. A real transfer function input is assumed to correspond to the magnitude response, and the associated phase response is estimated.

If FrequencyResponse is given, the alternative input fields ResonatorFrequency and ResonatorQL are ignored.

Optional, can be 'simulate' or 'compensate'. By default the effect of the resonator on the signal is simulated. If set to 'compensate' the pulse shape is adapted such that it compensates for the resonator profile by using a uniform critical adiabaticity criterion.

Relaxation can be controlled for the entire sequence or specific events:

Opt.Relaxation = true;   	% switches on relaxation for the entire simulation
Opt.Relaxation = [0 0 1];	% relaxation is active only during the third event

If relaxation is switched on, relaxation times (Sys.T1 and Sys.T2) have to be provided and the simulation is then performed in Liouville space, which slows down the simulation, for larger spin systems tremendously so. In order to avoid unexpected behavior, if the Opt.Relaxation is not used to switch relaxation on/off for the entire sequence, the length of it has to match the length of Exp.Sequence (all events need to be defined).

Optional, frequency of the (rotating) simulation frame in GHz. If SimFreq is not provided, a suitable simulation frame is automatically selected. No matter the choice of the simulation frame, all other frequencies must still be given in the lab frame - the program handles all required frequency shifts. For example, for running a Q-band simulation without using the default simulation frame, but one that rotates at 30 GHz:

Exp.mwFreq = 34.4;	        % Set carrier frequency to 34.4 GHz
Pulse.Frequency = [-50 50]; % Pulse with bandwidth of 100 MHz, sweeps from 34.35 to 34.45 GHz 

Opt.SimFreq = 30;  	      % Changes into a rotating frame at 30 GHz

To force a simulation in the labframe SimFreq has to be set to zero:

Opt.SimFreq = 0;  	      % Propagation in the labframe

Changing the simulation frame allows for using a larger time step, which in turn can strongly reduce computation time. Keep in mind, that if a time step is provided with Opt.IntTimeStep it is not adapted automatically. To have the program adapt the time step, remove the field Opt.IntTimeStep .

Integration time step in microseconds, optional. If no time step is provided, the program computes a suitable time step (taking into account the highest frequency that is being used for the pulse definition Exp.Frequency). Accurate results are obtained for an integration time step that corresponds to about 50 data points per oscillation or 1/50th of the largest time step size to fulfil the Nyquist criterion.

User provided integration time steps must at least fulfil the Nyquist criterion, or an error is thrown. For integration time steps larger than the recommended ones, a warning is printed to the command line.

If StateTrajectories is active, the output contains a cell array with density matrices at each time step during events that are detected. StateTrajectories can be switched on for the entire sequence as well as for specific events:

Opt.StateTrajectories = true;   	% switches on state trajectories for the entire simulation
Opt.StateTrajectories = [0 0 1];	% state trajectories is active only during the third event

State trajectories are only fully recorded for events where detection is active. Otherwise, the density matrices are stored only at the beginning and end of the event. In order to avoid unexpected behavior, if the Opt.StateTrajectories is not control the entire sequence, the length of it has to match the length of Exp.Sequence (all events need to be defined).

By default Ŝ_x is used as excitation operator. With the cell array ExcOperator it is possible to use custom excitation operators. The indexing of ExcOperator corresponds to the pulse index. It is possible to use the syntax from sop as well as matrices:

Opt.ExcOperator = {[] 'x(1|2)'};     % transition selective excitation operator for the second pulse
Opt.ExcOperator = {[0 1/2; 1/2 0]};  % custom excitation operator for the first pulse in matrix form

The imaginary part of a custom excitation operator is only considered if ComplexExcitation is active for the respective pulse.

Exp.nPoints = [3 4];
Exp.DetOperator = {'z1','+1'};

size(sig) = 

           3	  4	  11003		2 	

trace = sig(2,4,:,:);
% or if you want to remove the singular dimensions (e.g. for plotting)
trace = squeeze(sig(2,4,:,:));

size(t) = 

           1 	11003

size(t) = 

           3        4 	   11003

size(sig) = 

           3	  4

size(sig{1,1}) = 

           2 	11003
		   
size(sig{2,4}) = 

           2 	21003

size(t) = 

           3	  4
		   
size(t{1,1}) = 

           2 	11003
		   
size(t{2,4}) = 

           2 	21003

InterestingState = squeeze(out.FinalState(2,4,:,:)); % squeeze removes the singleton dimensions

InterestingState = 

		0.0592 - 0.0000i   0.4559 + 0.1965i
		0.4559 - 0.1965i  -0.0592 + 0.0000i

InterestingStateTrajectory = out.StateTrajectories{2,4}

size(out.StateTrajectories)

     3		4
size(InterestingStateTrajectory)

     1		21003