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


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;

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];


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];
Sys = nucspinadd(Sys,'63Cu',[50 500]);
A = [20 30];
Sys = nucspinadd(Sys,'14N',A);
Sys = nucspinadd(Sys,'14N',A);
Sys = nucspinadd(Sys,'14N',A);
Sys = nucspinadd(Sys,'14N',A);
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);
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];
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;

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

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

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

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).