Simulating solid-state cw EPR spectra

This user guide explains how to simulate solid-state cw EPR spectra of powders or single crystals using EasySpin's function pepper. It is assumed that you are familiar with the basics of MATLAB, esp. with structures.

It contains the following topics: There are the following advanced topics:
Running the simulation

Solid-state cw EPR spectra of powders and single crystals are computed by the EasySpin function pepper. It can be called with two or three parameters and can return both field axis and spectrum.

```pepper(Sys,Exp);                      % plots the spectrum
[Field,Spec] = pepper(Sys,Exp);       % returns the field axis and the spectrum
[Field,Spec] = pepper(Sys,Exp,Opt);   % additional simulation options in Opt
```

Don't forget the ; (semicolon) at the end of the line to suppress output to the screen.

The first argument `Sys` tells `pepper` all about the spin system, and the second argument `Exp` gives the experimental parameters. The third, optional, argument `Opt` contains settings concerning the simulation itself, like the number of orientations for powder simulations.

The outputs `Field` and `Spec` are arrays containing the magnetic field values and the spectrum, respectively. If no output is requested, `pepper` simple plots the spectrum. If the outputs are requested, `pepper` does not plot the spectrum, but you can plot it yourself using

```plot(Field,Spec);
```

Setting up a simulation and running it takes only a few lines of code. A very simple one would be

```Sys.g = [2 2.1];
Sys.lwpp = 0.5;
Exp.mwFreq = 9.5;
pepper(Sys,Exp);
```

This simulates and plots the 9.5 GHz EPR spectrum of an S=1/2 system with an axial g tensor. Copy and paste the code above to your MATLAB command window to see the graph.

The spin system

The first input argument to `pepper` is a structure specifying the spin system. It contains fields for the electron spin(s), the nuclear spins, and the various interaction parameters like g and hyperfine tensors.

`pepper` automatically assumes `S=1/2` for the spin quantum number. For systems with higher spin or more than one unpaired electron, the spin quantum number should be given in the field `S`.

```Sys.S = 1;          % a triplet state
Sys.S = 5/2;        % for e.g. high-spin Mn2+ or high-spin Fe3+
Sys.S = [1/2, 1/2];  % for a biradical
```

The field `g` contains the principal values of the g tensor(s). A simple rhombic S=1/2 system (e.g., a low-spin Fe3+) is

```Sys.g = [1.8, 2, 2.1];
```

Nuclear spins are included by specifying `Nucs` (comma-separated list of nuclei) and `A` (array of hyperfine tensor principal values, in MHz).

```Sys.Nucs = '2H';        % one 2H (deuterium) nucleus
Sys.A = [-1,-1,2]*4.2;  % hyperfine principal values in MHz
```

If the A tensor is tilted with respect to the molecular frame, the tilt angles can be provided via the field `AFrame`

```Sys.AFrame = [0 30 0]*pi/180; % [alpha beta gamma] in radians
```

The zero-field splitting is specified in the `D` field, in units of MHz. There are several different input possibilities:

```Sys.D = 120;          % D = 120 MHz, E = 0
Sys.D = [120 15];     % D = 120 MHz, E = 15 MHz
Sys.D = [-25,-55,80]; % principal values of D tensor, in MHz
```

D and E are related to the principal values of the D tensor (see reference page on the zero-field splitting).

Details about all the spin Hamilton parameters can be found on the spin Hamiltonian reference page. It is also possible to include several electron spins. Refer to the page about spin system structures for details.

No cw EPR spectrum is infinitely sharp. Lines are usually broadened due to several reasons. `pepper` provides means to include several line broadening models in a simulation.

The simplest way to include line broadening is to convolute a stick spectrum with a (Gaussian or Lorentzian) lineshape after the end of the simulation. Such a convolution broadening is specified in the spin system field `lwpp`.

```Sys.lwpp = 0.5;     % Gaussian broadening of 0.5 mT PP
Sys.lwpp = [0 2];   % Lorentzian broadening of 2 mT PP
Sys.lwpp = [1 2];   % Gaussian PP of 1mT + Lorentzian PP of 2 mT
```

The line width is in mT and refers to peak-to-peak (PP) line widhts. Instead, FWHM (full width at half height) line widths can be provided in the field `Sys.lw`.

```Sys.lw = 0.5;     % Gaussian broadening of 0.5 mT FWHM
Sys.lw = [0 2];   % Lorentzian broadening of 2 mT FWHM
```

For details about line shapes and conversion formulas to/from FWHM and peak-to-peak widths, see the page on line shapes.

Physically, there are several origins for line broadening. Large contribution to broadening often comes from unresolved hyperfine couplings and from distributions in the various magnetic parameters lige g, A and D that result from structural variations from one paramagnetic center to the next.

To include effects from unresolved hyperfine couplings, an orientation-dependent phenomenological broadening can be specified in `HStrain`:

```Sys.HStrain = [50 50 87];   % [along x, along y, along z], in MHz
```

Distributions in magnetic parameters are called strains. g and A strains are given in similar fields:

```Sys.gStrain = [0.01 0.02 0.005];
Sys.AStrain = [10 10 30]; % in MHz
```

The three values in `gStrain` are the FWHM parameters of the Gaussian distributions of the respective g principal values given in `Sys.g`. `AStrain` is the same for the A tensor. The g and A strains are correlated.

Distributions of the D tensor values can be given in `DStrain`, where the first value is the width of the (scalar) D distribution, and the second is the width for the E distribution.

All these broadening parameters can be combined. However, usually a modelling of the broadening with `lwpp` or `HStrain` is absolutely sufficient.

Basic experimental settings

The second input argument, `Exp`, collects all experimental settings. Just as the spin system, `Exp` is a structure containing several fields.

Microwave frequency. To simulate an EPR spectrum, Easyspin needs at a minimum the spectrometer frequency. Put it into `Exp.mwFreq`, in units of GHz.

```Exp.mwFreq = 9.385;  % X-band
Exp.mwFreq = 34.9;   % Q-band
```

Field range. There are two ways to enter the magnetic field sweep range. Either give the center field and the sweep width (in mT) in `Exp.CenterSweep`, or specify the lower and upper limit of the sweep range (again in mT) in `Exp.Range`.

```Exp.CenterSweep = [340 80]; % in mT
Exp.Range = [300 380];      % in mT
```

On many cw EPR spectrometers, the field range is specified using center field and sweep width, so `Exp.CenterSweep` is often the more natural choice.

`Exp.CenterSweep` and `Exp.Range` are only optional. If both are omitted, EasySpin tries to determine a field range large enough to accomodate the full spectrum. This automatic ranging works for most common systems, but fails in some complicated situations. EasySpin will issue an error when it fails.

Points. By default, `pepper` computes a 1024-point spectrum. However, you can change the number of points to a different value using

```Exp.nPoints = 5001;
```

You can set any value, unlike some EPR spectrometers, where often only powers of 2 are available (1024, 2048, 4096, 8192).

Harmonic. By default, EasySpin computes the first-harmonic absorption spectrum, i.e. the first derivative of the absorption spectrum. By changing `Exp.Harmonic`, you can request the absorption spectrum directly or the second-harmonic (second derivative) of it.

```Exp.Harmonic = 0; % absorption spectrum, direct detection
Exp.Harmonic = 1; % first harmonic (default)
Exp.Harmonic = 2; % second harmonic
```

Modulation amplitude. If you want to include effects of field modulation like overmodulation, use `Exp.ModAmp`

```Exp.ModAmp = 0.2; % 0.2 mT (2 G) modulation amplitude, peak-to-peak
```

Time constant. To include the effect of the time constant, apply the function rcfilt to the simulated spectrum.

More experimental settings

For more advanced spectral simulations, EasySpin offers more possibilities in the experimental parameter structure `Exp`.

Mode. Most cw EPR resonators operate in perpendicular mode, i.e., the oscillating magnetic field component of the microwave in the resonator is perpendicular to the static field. Some resonators can operate in parallel mode, where the microwave field is parallel to the static one. EasySpin can simulate both types of spectra:

```Exp.Mode = 'perpendicular'; % perpendicular mode (default)
Exp.Mode = 'parallel';      % parallel mode
```

Temperature. The polarizing effects of low sample temperatures can also be included in the simulation by specifying the temperature:

```Exp.Temperature = 4.2; % temperature in kelvin
```

With this setting, EasySpin will include the relevant polarization factors resulting from a thermal equilibrium population of the energy levels. For S=1/2 systems, it is not necessary to include the temperature. However, it is important in high-spin systems with large zero-field splittings, and in coupled spin systems with exchange couplings.

Microwave phase. Occasionally, the EPR absorption signal has a small admixture of the dispersion signal. This happens for example when the microwave phase in the reference arm is not absolutely correctly adjusted. EasySpin can mix dispersion with absorption if a Lorentzian broadening is given:

```Sys.lwpp = [0.2 0.01];           % Lorentzian broadening (2nd number) required

Exp.mwPhase = 0;                 % pure absorption
Exp.mwPhase = pi/2;              % pure dispersion
Exp.mwPhase = 3*pi/180;          % 3 degrees dispersion admixed to absorption
```
Powders and crystals

If not specified otherwise, `pepper` computes a powder spectrum. But it is as well straightforward to simulate spectra for a single crystal. The orientation (or orientations if more than one) of the single crystal can be provided in the experiment structure field `Exp.CrystalOrientation`. This field should contain the tilt angles between crystal and laboratory frame (right-handed coordinate system with z along the static field and x along the microwave magnetic field), one set of three angles per row.

For a crystal with its molecular frame aligned with the laboratory frame, the setting is

```Exp.CrystalOrientation = [0 0 0];
```

If you need more than one crystal at the same time, then just specify more than one orientation.

```Exp.CrystalOrientation = [0 0 0; 0 pi/4 0];
```

If `Exp.CrystalOrientation` is missing or set to `[]` (an empty array), `pepper` simulates the powder spectrum.

In many crystals, there are several symmetry-related sites with identical paramagnetic centers differing only in their orientation in the crystal. You can tell `pepper` about this by providing the crystal symmetry in the field `Exp.CrystalSymmetry`, e.g.

```Exp.CrystalSymmetry = 'P21/m'; % space group symbol
Exp.CrystalSymmetry = 11;      % space group number (between 1 and 230)
Exp.CrystalSymmetry = 'C2h';   % point group, Schönflies notation
Exp.CrystalSymmetry = '2/m';   % point group, Hermann-Mauguin notation
```

With the crystal symmetry given, `pepper` not only computes the spectrum for the orientation given in `Exp.CrystalOrientation`, but also for all symmetry-related sites.

Simulation options

The third input argument to `pepper` contains simulation options. All of them have reasonable default values, but sometimes it might be necessary to change one of them. In the following the most important ones are presented.

If you want `pepper` to print information about the simulation to the command window during the computation, use

```Options.Verbosity = 1;
```

`'Verbosity'` tells `pepper` how much of progress information to show in the command window. 0 (the default) suppresses all output, 1 is normal intormation, and 2 prints more information, relevant only for debugging.

Another useful option is `nKnots`, which determines how many orientations `pepper` will include in the simulation of a powder spectrum. If this value is too low, the spectrum shape contains ripples. `nKnots` is the number of orientations between the z axis and the x axis (between theta = 0 and theta = 90 degrees).

```Options.nKnots = 31; % corresponds to 3-degree increments
```

The higher `nKnots`, the finer the orientational grid. The default value of 19 (5-degree increments) is appropriate for most systems. A value larger than 181 (0.5-degree increments) is rarely needed.

After having computed the spectrum for a number of orientations specified by `nKnots`, the simulation function interpolates these spectra for additional orientations before summing up all spectra. This interpolative refinement can be configured with a second number in `nKnots`. `nKnots = [19 4]` means that `pepper` interpolates additional 4 spectra between two adjacent orientations evaluated.

```Options.nKnots = [19 10];  % massive interpolation
Options.nKnots = [19 0];   % no interpolation
```

The option `Output` can be used to determine the form in which `pepper` returns the spectral data.

```% single crystal: orientations separately
% powders: transitions separately
Options.Output = 'separate';

% sum over all orientations and transitions
Options.Output = 'summed';
```

There are more option fields available. For details, see the documentation page on pepper.

Multiple components

Often, an EPR spectrum shows a mixture of spin species. To simulate these spectra, each of the component spectra has to be simulated and added with the appropriate weight (depending on spin concentration) to the total spectrum.

This can be done automatically by `pepper`. Just provide the component spin systems with their weights as a cell array (in curly braces) to `pepper`. For example, here is the simulation of a very simple two-component mixture with 2:1 ratio of spin concentrations.

```Sys1.g = 2;
Sys1.lwpp = 1;
Sys1.weight = 2;
Sys2.g = 2.1;
Sys2.lwpp = 0.8;
Sys2.weight = 1;

Exp.mwFreq = 9.5;
Exp.Range = [300 360];

pepper({Sys1,Sys2},Exp);
```

You don't have to specify `Sys.weight` - if it's absent it is assumed to be 1. These weights are absolute, i.e. a simulation with `Sys.weight=20` yields a spectrum that is 10 times more intense than the one obtained with `Sys.weight=2`. There is no limit to the number of components in a simulation.

Systems with several nuclei

`pepper` uses matrix diagonalization as thed default method for simulating spectra. For systems with several nuclei this can be very time-consuming. To accelerate such computations, `pepper` provides first- and second-order perturbation theory as an alternative methods. The relevant simulation option that tells EasySpin about is `Opt.Method`.

As an example, we look at the simulation of the spectrum of Cu2+ porphyrin.

```Sys.S = 1/2;
Sys.g = [2 2.2];
A = [20 30];
Sys.lwpp = 0.5;
```

With matrix diagonalization (`Opt.Method='matrix'`, which is the default), the simulation needs several hours. With second-order perturbation theory (`Opt.Method='perturb2'`), the simulation is orders of magnitude faster, but potentially less accurate. We can compare the full matrix diagonalization to the perturbation simulation.

```Exp.mwFreq = 9.5;
Exp.Range = [260 380];
Opt.Method = 'matrix';
[x,y1] = pepper(Sys,Exp,Opt);
Opt.Method = 'perturb2';
[x,y2] = pepper(Sys,Exp,Opt);
plot(x,y1,'r',x,y2,'b');
```
Non-equilibrium populations

`pepper` can handle both thermal equilibrium and non-equilibrium populations. Both are specified in the field `Temperature` of the experimental settings structure.

For thermal equilibrium, just give the temperature in kelvin:

```Exp.Temperature = 77; % 77K, boiling point of liquid nitrogen
```

For non-equilibirum populations, `Temperature` must be a vector. If the spin systems contains N electron states, then this vector must contain N elements, each specifying the population of one of the electron states at zero field, sorted according to their energy from lowest to highest.

E.g., an organic triplet with S=1 and I=1 has 3 electron states, each further split into three sublevels by the coupling to the nuclear spin. The population vector in this case should contain three elements:

```Exp.Temperature = [0.6 0.8 1.1]; % highest state is most populated
```

This specifies that all the sublevels of the lowest zero-field electron states have a population of 0.6, etc. The sublevels of the highest-energy zero-field electron state have a population of 1.1. The populations don't have to be normalized, `pepper` takes care about that.

To compute the state populations for a non-zero field state, `pepper` decomposes it into a linear combination of zero-field states and combines the zero-field populations using the resulting linear combination coefficients.

A simple example of a non-equilibrium triplet system is

```Sys.S = 1; Sys.g = 2; Sys.lw = 0.2;
Sys.D = 100;
Exp.mwFreq = 9.5; Exp.Range = [320 360]; Exp.Harmonic = 0;
Exp.Temperature = [0.5 0.6 0.9];
pepper(Sys,Exp);
```
Partially ordered systems

In powders and frozen solutions (disordered systems), paramagnetic molecules are randomly oriented in the sense that any orientation can occur with equal probability. In other situations, like in polymers, biomembranes or liquid crystals, the paramagnetic molecules may be partially aligned or ordered, so that some orientations are more probable than others. As a result, the spectra of such partially ordered systems are different from those of powders and frozen solutions.

pepper can include partial ordering in the spectral simulation. For this, set a value (different from zero) in the experiment structure field `Exp.Ordering`. It is a number which specifies the nature and the degree of the ordering.

If it is zero, no ordering is used. Positive values mean that the molecules are partially aligned so that the magnetic field vector is mostly approximately parallel to the molecular z axis. Negative values specify partial ordering such that molecular orientations with the magnetic field vector in the molecular xy plane are more probable.

```Exp.Ordering = 0;   % all orientations equally populated
Exp.Ordering = -1;  % slightly preferential orientation with field in the molecular xy plane
Exp.Ordering = +10; % strongly aligned such that field is along the molecular z axis
```

The ordering function used is very simple (see the pepper documentation for more information), but sufficient for many cases. Here is a simulation of a sample where the molecules are preferentially oriented such that the molecular z axis is close to the magnetic field vector:

```Sys.g = [2 2 2.2];
Sys.lwpp = 1;
Exp.mwFreq = 9.5;
Exp.Ordering = +2;
pepper(Sys,Exp);
```

You can also define your own custom orientational distribution in a separate function and supply to `pepper` as a function handle in `Exp.Ordering`. See the pepper documentation for details. Written as an anonymous function, the built-in orientational distribution is equivalent to

```Exp.Ordering = @(phi,theta) exp(lambda*plegendre(2,0,cos(theta)));
```

where the EasySpin function plegendre returns the associated Legendre polynomial and `lambda` corresponds to the number given in `Exp.Ordering`.

Frequency-swept spectra

`pepper`, like the other cw EPR simulation functions `garlic` and `chili`, does field sweeps by default. However, you can it them to simulate frequency-swept spectra as well.

For this, all you need to do is the following

• Give a static magnetic field (in mT) in `Exp.Field`. Make sure you do not set `Exp.mwFreq`, otherwise EasySpin does not know what to do.
• Give a frequency range (in GHz) in `Exp.mwRange` or `Exp.CenterSweep`. You can also omit these, in which case `pepper` will determine an adequate range automatically.
• If you use `Sys.lw` or `Sys.lwpp`, make sure they are in MHz units. For a frequency sweep, these convolutional linewidth parameters are understood to be in MHz (and not in mT, as they are for field sweeps).

Here is an example of a frequency-swept spectrum of a nitroxide radical:

```clear
Sys.g = [2.008 2.006 2.002];
Sys.Nucs = '14N';
Sys.A = [20 20 100];
Sys.lwpp = 10;           % peak-to-peak linewidth, in MHz
Exp.Field = 340;         % static field, in mT
Exp.mwRange = [9.3 9.7];   % frequency range, in GHz
pepper(Sys,Exp);
```

By default, `pepper` returns the absorption spectrum (`Exp.Harmonic=0`) when you simulate a frequency-swept spectrum. To get the first or second derivative, change `Exp.Harmonic` to 1 or 2. Note however that `Exp.ModAmp` is not supported for frequency sweeps.

All other capabilities of `pepper` apply equally to frequency sweep and to field sweeps. For example, you can simulate crystals or multi-component spectra, and you can change the excitation mode. Importantly, just like field sweeps, frequency sweeps can be simulated using different methods: matrix diagonalization (`Opt.Method='matrix'`; highly accurate but potentially slow) or perturbation theory (`Opt.Method='perturb'`, `Opt.Method='perturb1'`, `Opt.Method='perturb2'`; less accurate but much faster).