Simulation of spin dynamics with arbitrary wave forms
spidyan(Sys,Exp) spidyan(Sys,Exp,Opt) sig = spidyan(...) [t,sig] = spidyan(...) [t,sig,out] = spidyan(...)
See also the examples on how to use spidyan
and the userguide.
This function simulates spin dynamics of pulse EPR experiments with arbitrary waveforms.
The focus of spidyan
lies on investigation of spin dynamics and on effects of non-ideal or frequency-swept pulses on the outcome of experiments.
Hence, the outputs of spidyan
are a time domain signal and density matrices of a single spin packet.
For simulating transients or spectra of realistic samples, use the function saffron
, which is based one the same propagation engine.
Time-domain signals or state trajectories simulated by spidyan
are typically further analyzed by (Matlab) software home-written by the user.
At the expense of requiring more knowledge on spin quantum mechanics than higher-level EasySpin functions, spidyan
provides a lot of freedom in specifying input and output (e.g., custom (non-physical) initial states, transition selective excitation operators, freely-choosable detection operators,...).
This freedom should help you to investigate your spin dynamics problem in depth.
There are up to three possible output arguments.
t
(time in microseconds) and the simulated data in sig
(time domain trace).
For a single acquisition or a multidimensional experiment where all time traces have the same length, t
and sig
are numeric arrays.
For multidimensional simulations where the length of the trace changes, t
and sig
are cell arrays.
out
contains additional output such as the final state in out.FinalState
(or states for a multidimensional experiment) after propagation and the cell array out.StateTrajectories
which contains state trajectories (density matrices) at each point of detection, but is available only if requested through Opt.StateTrajectories
.
A detailed discussion of the structure of the output can be found further below.
The three input arguments to spidyan
are
Sys
: spin system (paramagnetic molecule)Exp
: experimental parametersOpt
: simulation options
Sys
is a spin system structure.
Most of the regular fields of Sys
can be used for the construction of the spin Hamiltonian.
Line broadening parameters used by other simulation functions (lw
, lwpp
, gStrain
, etc.) are ignored.
Some additional fields can be used with spidyan
: To facilitate working with frequency-swept pulses, the field ZeemanFreq
is available.
For simulations with relaxation, at least generic relaxation times T1
or T2
should be provided.
If only one is provided, the other one is set by default to a very large value (1010 microseconds), essentially switching that type of relaxation off.
ZeemanFreq
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.
T1
spidyan
to simulate relaxation effects, it is necessary to switch relaxation on with Opt.Relaxation
(see below).
T2
spidyan
to simulate relaxation effects, it is necessary to switch relaxation on with Opt.Relaxation
(see below).
initState
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
eqState
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.
Exp
contains experimental parameters.
Field
Sys.g
.
If Sys.ZeemanFreq
is used, Field
can be omitted.
Sequence
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 microsecondsThe delays are defined to go from the end of one previous element in
Sequence
to the beginning of the next (unlike in Bruker spectrometers).
mwFreq
Pulse.Frequency
) need to be defined relative to it:
Exp.mwFreq = 33.5; Pulse.Frequency = [-250 250]; % 500 MHz sweep
PhaseCycle
PC = [0, 1; pi, -1]; % phase cycle Exp.PhaseCycle = {[] PC}; % phase cycles the second pulsePhase 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
).
mwPolarization
'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.
Detection:
DetWindow
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 delayIn 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 tauOf 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 nsIf 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
DetSequence
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 detectionBesides 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 detectedIf 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.
DetOperator
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 spinDetection 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 formFor 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.
DetFreq
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 shiftedFor 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
DetPhase
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.
By defining nPoints
and DimX
it is possible to vary parameters of the pulse sequence and create one or multidimensional experiments.
This also changes the structure of the output compared to a single acquisition.
More details on how the output changes can be found below.
nPoints
Exp.nPoints = [10 150]; % 10 points along the first indirect dimension, 150 along the second
Dim1
, Dim2
,... 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 stepFor 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 stepSeveral 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 steppedFor 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 dimensionA 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 delayHowever, 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 stepAn 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 stepThis 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 sweepA 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 stepIn 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
Bandwidth limiting effects of a resonator on a pulse can be incorporated in two ways:
The resonator profile is defined using the resonator center frequency ResonatorFrequency
and the loaded Q-value ResonatorQL
.
The resonator frequency response is then computed from the ideal transfer function for an RLC series circuit)
Exp.ResonatorFrequency = 33.5; % resonator center frequency in GHz Exp.ResonatorQL = 300; % loaded Q-value
Or the transfer function in combination with the frequency axis are given with FrequencyResponse
.
Exp.FrequencyResponse = [Frequency; TransferFunction] % transfer function and its correspond frequency axis
By setting ResonatorMode
it is possible to compensate for the resonator.
Compensation aims to provide a waveform that excites spin packets in the resonator with constant critical adiabaticity if the originally specified waveform would lead to excitation with constant critical adiabaticity in the absence of a resonator (mainly for chirp and hyperbolic secant pulses).
For more details see pulse and resonator.
ResonatorFrequency
ResonatorQL
FrequencyResponse
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.
ResonatorMode
'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.
Opt
contains additional simulation parameters.
Relaxation
Opt.Relaxation = true; % switches on relaxation for the entire simulation Opt.Relaxation = [0 0 1]; % relaxation is active only during the third eventIf 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).
SimFreq
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 GHzTo force a simulation in the labframe
SimFreq
has to be set to zero:
Opt.SimFreq = 0; % Propagation in the labframeChanging 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
.
IntTimeStep
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.
StateTrajectories
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 eventState 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).
ExcOperator
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 formThe imaginary part of a custom excitation operator is only considered if
ComplexExcitation
is active for the respective pulse.
The data type of the returned time-domain signal depends on experiment parameters.
If you run a single experiment (no Exp.nPoints
), the output sig
is a two-dimensional numeric array, with the first dimension corresponding being the detected transient and the second dimension specified detection operators.
If one the specified detection operators is a ladder operator, say Ŝ-, sig
is a complex-valued vector.
The time axis is returned in t
and will be a vector with the same length as sig
.
If you run a multidimensional experiment, sig
can be an (n+2)-dimensional numeric array (if each acquisition point has the same length) or an n-dimensional cell array (if the detection length changes), where n is the number of indirect dimensions of your simulation.
The indexing of n corresponds to how Exp.nPoints
was defined:
If the total length of all detected events is identical for each acquisition point, sig
will be a (n+2) dimensional array.
For Exp.nPoints = [3 4];
a total of 12 data sets will be returned in a four dimensional array:
Exp.nPoints = [3 4]; Exp.DetOperator = {'z1','+1'};
size(sig) = 3 4 11003 2The first dimension of your output has a length of
3
and the second dimension 4
.
These correspond to the indirect dimensions.
The second-to-last dimension corresponds are the expectation values at each point of time during which was detected (the transient) and the last dimension to the detection operators.
To get the time trace that corresponds to the second point in the first dimension and the last point in the second dimension you could use:
trace = sig(2,4,:,:); % or if you want to remove the singular dimensions (e.g. for plotting) trace = squeeze(sig(2,4,:,:));If all acquisition points have the same time axis,
t
is a vector:
size(t) = 1 11003If not,
t
is an (n+1) dimensional numeric array, with the same indexing for n as for sig
.
This can be encountered when detecting a moving echo, with a detection window that is always centered around the echo, but moves in terms of absolute time (see example on two-pulse ESEEM).
size(t) = 3 4 11003
When the detection length changes (e.g. if you vary the length of a pulse or a free evolution event that is being detected) between acquisition points, t
and sig
become n-dimensional cell arrays.
The indexing of the elements in each cell array again corresponds to Exp.nPoints
:
size(sig) = 3 4Each element of
sig
is a two-dimensional numeric array with the first dimension corresponding to the detection operators and the second to the time axis:
size(sig{1,1}) = 2 11003 size(sig{2,4}) = 2 21003The time axes for each acquisition point are stored as vectors in the elements of
t
:
size(t) = 3 4 size(t{1,1}) = 2 11003 size(t{2,4}) = 2 21003
For a single acquisition, the output out.FinalState
is the final density matrix of size nH x nH, where nH is the dimensionality of your spin system in Hilbert space.
If a simulation with more than one acquisition point is run, s
is an (n+2)-dimensional numeric array, with the same rules for indexing as demonstrated above for sig
.
E. g. to get the final state of the second data point in the first dimension and the fourth data point in the second dimension, write:
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
If Opt.StateTrajectories
is active during at least one event, out.StateTrajectories
contains the state trajectories.
This means, it contains the state (density matrix) of each propagation step during the events selected in Opt.StateTrajectories
.
out.StateTrajectories
is an n-dimensional cell array where the dimensions once again correspond to Exp.nPoints
.
To get the state trajectories for the second data point in the first dimension and for the fourth datapoint in the second dimension of the example above, write:
InterestingStateTrajectory = out.StateTrajectories{2,4}
size(out.StateTrajectories) 3 4 size(InterestingStateTrajectory) 1 21003For full recording of state trajectories, detection needs to be active during those events. Otherwise only the states at the beginning and the end of that event are recorded.
saffron, resonator, signalprocessing, pulse, sop