Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Appl Magn Reson. Author manuscript; available in PMC 2010 December 1.
Published in final edited form as:
Appl Magn Reson. 2009 December 1; 36(2-4): 237–258.
doi:  10.1007/s00723-009-0023-5
PMCID: PMC2786184

Calculation of Double-Quantum-Coherence Two-dimensional Spectra: Distance Measurements and Orientational Correlations


The double quantum coherence (DQC) echo signal for two coupled nitroxides separated by distances 10 Å, is calculated rigorously for the six-pulse sequence. Successive application of six pulses on the initial density matrix, with appropriate inter-pulse time evolution and coherence pathway selection leaves only the coherent pathways of interest. The amplitude of the echo signal following the last π pulse can be used to obtain a one-dimensional dipolar spectrum (Pake doublet), and the echo envelope can be used to construct the two-dimensional DQC spectrum. The calculations are carried out using the product space spanned by the two electron-spin magnetic quantum numbers m1, m2 and the two nuclear-spin magnetic quantum numbers M1, M2, describing e.g. two coupled nitroxides in bilabeled proteins. The density matrix is subjected to a cascade of unitary transformations taking into account dipolar and electron exchange interactions during each pulse and during the evolution in the absence of a pulse. The unitary transformations use the eigensystem of the effective spin-Hamiltonians obtained by numerical matrix diagonalization. Simulations are carried out for a range of dipolar interactions, D, and microwave magnetic field strength B for both fixed and random orientations of the two 14N (and 15N) nitroxides. Relaxation effects were not included. Several examples of one- and two-dimensional Fourier transforms of the time domain signals vs. dipolar evolution and spin-echo envelope time variables are shown for illustration. Comparisons are made between 1D rigorous simulations and analytical approximations. The rigorous simulations presented here provide insights into DQC ESR spectroscopy, they serve as a standard to evaluate the results of approximate theories, and they can be employed to plan future DQC experiments.

1. Introduction

The first theoretical analysis of DQC spectra was undertaken by Saxena and Freed (SF) [1]. They briefly summarized results for constant-time 6-pulse DQC sequence, but focused instead on a “forbidden” DQC 5-pulse sequence. In fact, it is the 6-pulse sequence, and not the 5-pulse sequence of that type, which has proved to be most useful in distance measurements [2-5]. While the equations in [1] are generally correct, their prediction that the “allowed” DQC signals would be very weak, have been contradicted by the strong DQC signals obtained in many experimental studies and their analysis [2-5]. It is presumed that there were undetected problems with the numerical simulations.

When a sample containing bilabeled proteins is subjected to a sufficiently strong microwave pulse, the nitroxide ESR spectrum is almost uniformly excited, so that any orientational selection is largely suppressed, that is, it does not modify the echo amplitude (except for the effect of pseudosecular dipolar terms, essential for short distances). Also, as we show, in high B1-fields (B12D), the effect of dipolar coupling during the pulses becomes relatively weak. Therefore, for not very short distances and in sufficiently strong B1's, the information on orientations of the magnetic tensors of the spin-label moieties, is virtually excluded from the time-domain dipolar evolution of the echo amplitude, taken at its maximum. However, as we show, it is still retained in the spin-echo envelope and can be retrieved by recording the 2D time-domain data as a function of the spin-echo time (techo) and the dipolar evolution time (tdip) and then converting it into a 2D-FT spectrum, which, after making a “shear” transformation [6] separates the dipolar dimension from the spectral dimension. Rigorous computations of 1D and 2D signals have been carried out and are presented here. Efficient but approximate analytical expressions to this end were developed for 1D signals by Borbat and Freed [2] (cf. Appendix III), who omitted the dipolar coupling during the pulses and assumed an ideal double-quantum (DQ) filter. Whereas such expressions are quite useful for practical purposes and computationally very efficient, these approximations may not be generally valid, especially in the case of short distances (e.g. less than 15.0 Å). In order to test the nature and extent of deviations from the exact results and to establish the scope of applicability, numerical simulations of 1D spectra were carried out rigorously using the new codes developed for 2D computations. Unlike [1], the pulse propagators are calculated, using highly accurate numerical diagonalizations of the Hamiltonians involved, rather than applying a Trotter expansion [7, 8]. Although the computational approach is necessarily time-consuming, it does provide deeper insights into the features of DQC spectroscopy.

2. Theoretical Background

The pulse pattern and the relevant coherence pathways for the 6-pulse DQC sequence are shown in Fig. 1 [2, 3]. The method used to compute the 2D-DQC spectra is outlined as follows.

Figure 1
In this diagram (a) shows the 6-pulse DQC sequence. The coherence pathways in (b) correspond with the pulses shown in (a) in that a transition from one p state to another p state is generated by a pulse; the horizontal lines show coherence orders during ...

The initial density matrix operator in thermal equilibrium for the two nitroxides described by the static spin-Hamiltonian Ĥ is:


where the z-axis is defined to be aligned along the direction of the external magnetic field, and the subscripts number the two electron spins. The time evolution of the spin density matrix, ρ(t), is governed by the Liouville-von Neumann equation [1]:


where H^^ρ[H^,ρ], ħ is Plank's constant divided by 2π, Γ^^ is the relaxation operator, i2 = −1. Neglecting the relaxation, the density matrix evolves, under the action of Ĥ as:


Thus after a period of time Δt the density matrix, ρ(t) becomes:


A numerical implementation of Eq. 4 to compute the DQC ESR signal is outlined below. Note that in the sequel the Hamiltonians will be expressed in angular frequency units.

3. Computation of echo signal

One starts with the (un-normalized) initial density matrix in thermal equilibrium:


for the two coupled nitroxides, with electron spins Sk = 1/2, (normalization will be performed at the end of the calculation). The 6-pulse sequence [9], as shown in Fig. 1, transforms the initial density matrix under the successive action of six pulse propagators Rk (k=1, 2, ... 6), due to the pulses, and six subsequent free-evolution propagators Qk (k = 1, 2,... ,6) . This sequence can be defined as follows:


The twelve time evolution periods described by Eq. 6 lead to the density matrix in the final form ρf. This is calculated as follows. A k-th pulse, applied at the time t and acting during the period of time, τk, in the frame of reference rotating with the angular frequency of the circular component of microwave magnetic field resonant with Larmor frequency of the nitroxide electron spin, transforms the density matrix, ρ(t), according to:


or in the notation of Eq. 6,


with Rk is the k-th pulse propagator due to the effective Hamiltonian Hk acting during the period of time τk. The action of a π-pulse can change the sign of a coherence order, p, and the π/2 pulse can generate other coherence orders [10]. In order to follow the pathways of interest, the density matrix is then projected onto the coherence pathways pk, which are chosen after the pulse according to Fig. 1, as follows:


where the idempotent operator P (pk) projects the density matrix on the coherence pathways pk chosen after the k-th pulse. As shown in Fig. 1, the coherence pathways of interest are [(1, −1); (−1, 1); (2, −2); (−2, 2); (1); (−1)] after the actions of the six pulses, with two branching points leading to a total of four distinct pathways. The subsequent free evolution during the time tk transforms ρ′(t + τk) as ρ(t+τk)Q(tk)ρ(t+τk+tk) according to:


where H is the spin-Hamiltonian in the absence of a pulse. The density matrix ρ(t + τk + tk) is then used in place of ρ(t) in Eq. 9 for the calculation of the density matrix under the action of the (k+1)-pulse, and the steps defined by Eqs. 910 are repeated to arrive at the final density matrix ρfρ(k(τk+tk)). Computed in this way, the final density matrix thus becomes a function of several arguments, ρf = ρ(τ, t, p, tdip, techo, η, λ1, λ2) . The arguments are as follows: τ = (τ1,...,τ6) are the pulse durations; t = (t1,...,t6) are the subsequent free evolution periods; p = (p1,...,p6) are the relevant coherence orders during the evolution periods; tdip, techo are time variables used to record the dipolar evolution and to produce the spin echo envelope. The remaining arguments are the Euler angles η = (χ, θ, [var phi]), which define in the laboratory frame the orientation of the vector r connecting the magnetic dipoles associated with the electron spins; and in this dipolar (molecular) frame, whose z-axis is coincident with r, the Euler angles λk = (αk, βk, γk) define the principal axis of the nitroxide magnetic tensors (Fig. 2) with α1 chosen to be 0. (The angle χ was set to zero as is appropriate for isotropic media). Finally, the complex echo signal is given by:


where ρ~f is the normalized density matrix.

Figure 2
The set of Euler angles λk = (αk, βk, γk), (k = 1, 2), which define the orientations of the hf and g-tensor tensor principal axes for nitroxides 1 and 2 in the dipolar (molecular) frame of reference. In this frame the ...

The evolution thus depends on the exact form of various propagators that is given by H0 in the absence of a pulse or H0 + Hp in the presence of a pulse. Appropriate pulse time intervals τk are chosen to achieve nominal flip angles of π/2 ( k=1,3,5) and π(k=2,4,6), respectively. Here,




for k=1,2 denoting nitroxides 1 and 2 and H12 giving their coupling


Here, J is the electron exchange constant and D is the dipolar coupling constant


where γe is the gyromagnetic ratio for an electron, h is Plank's constant, and r is the distance between the nitroxide's magnetic dipoles separated by r and considered in the point dipole approximation. (The dipolar constant d = 2D/3 will most often be used throughout the text). In Eq. 13, I1,2 are the nuclear spins of the nitrogen (14N or 15N) nuclei on the two nitroxides.

The interaction with the radiation field due to the applied microwave pulse k in the reference frame rotating with the carrier frequency ωrf (that is usually set at or near the Larmor frequency) is given by:


where B1k is the amplitude of the circular magnetic component. The phases, [var phi]k, can be set to zero for all the pulses for purposes of the present calculations and consequently H pk = γeB1kSx. The amplitudes, B1k, do not have to be equal for different k's, but we will show results for the simplest case of equal amplitudes. H0k (cf .Appendix I) can be written in the laboratory frame as


The two-dimensional time-domain DQC signal is calculated for a given set of λk and η, using appropriate variations of the time intervals following the various pulses; as given in the caption of Fig. 1 for typical values used for simulation. The calculations were carried out in the product space S1 [multiply sign in circle] S2 [multiply sign in circle] I1 [multiply sign in circle] I2 with the dimension N=(2S1+1)(2S2+1)(2I1+1)(2I2+1) (or using those generated by subspace permutations, as needed). H0 and H are represented by order N×N matrices H0 and H, i.e. with the size of 36×36 for the two coupled (14N) nitroxides, (S1,2=1/2, I1,2=1). H0, not including H12, can be conveniently diagonalized in the product space S1 [multiply sign in circle] I1 [multiply sign in circle] S2 [multiply sign in circle] I2 by the unitary transformation


wherein S, S are the matrices representing polarization operators for spin 1/2 and T, T are dimension (2Ik+1) unitary matrices diagonalizing the respective nuclear manifolds, α and β. The transformation, however, brings Skx = Sx [multiply sign in circle] 1Ik to the form, which contains off-diagonal elements in transformed 1Ik's.


where Mk=TkαTkβ. If nuclear Zeeman terms can be neglected, as is appropriate for nitroxides in magnetic fields of up to about 12 kG, which corresponds to Q band, Mk=1Ik. In this case both H1 and H12 become block-diagonal in I1 [multiply sign in circle] I2 [multiply sign in circle] S1 [multiply sign in circle] S2 with (2I1+1)(2I2+1) blocks, corresponding to order 4×4 S1 [multiply sign in circle] S2 electronic subspace. This block-diagonal form makes computations run significantly faster. At higher fields nuclear Zeeman terms should be retained, leading to coupling between nuclear and electron coherences and a large increase in computation time. The procedure to calculate ρf(tdip,techo,θ,[var phi]) is outlined in Appendix II.

Finally, the complex echo signal is given by:


where 1N is the unity matrix in the product space. The echo signal from a powder (or macroscopically-aligned) sample is the average of the signals over the orientations of the molecule in the laboratory frame:


Here PΩ is the angular distribution of molecular axes in the laboratory frame (PΩ=1/4π for an isotropic distribution). In performing powder averaging in isotropic medium it suffices to set the integration limits to [0, π] in axial angles and [0, π/2] in polar angles. (Anisotropic media will require averaging over all three Euler angles in η). When needed, averaging over some or all of the remaining Euler angles, λk, and over distances can be conducted.

4. Computational Efficiency and Features

The 1D and 2D simulations based on block-diagonal approximations (i.e. after neglecting the nuclear Zeeman terms in H0) were carried out with a typical PC using a software code written in MATLAB. In a magnetic field of 6,200 G (corresponding to Ku band) the outcome was virtually indistinguishable from that produced by the most rigorous simulations, and the examples we show were produced in this way. Powder averaging employed averaging over the grid on unit sphere (actually over an octant), but also was carried out based on Monte-Carlo averaging. 2D simulations used typically 400×80 points in the time-domain, a 180×90 point (or greater) grid in {cosθ, [var phi]} or {θ, [var phi]}, or else 1−4·104 Monte-Carlo trials. For 1D simulations, the grids can be significantly larger. The grid in θ (or cosθ) was spiral or else randomized within the Δθ (or Δcosθ) limits on every step in [var phi]. Grids were constructed for the cosines or the angles for the other polar angles, βk. The MATLAB code is general enough to permit averaging over any subset of Euler angles or indeed over all of them whether in mesh mode or in Monte Carlo averaging. Averaging was also carried out with respect to the dipolar constant using suitable distance distributions, P(r). The Monte Carlo mode also included angular-dependent inhomogeneous broadening. The last two features are the most practical for 1D computation. A typical computation time for 2D simulations is about 1−2 hours, but could be as long as 10 hours for very precise results that require the largest meshes noted above or for long Monte Carlo averaging, (e.g. when d is large or distribution in d is broad). A very precise 1D simulation requires some 10−50 min, whereas 1D simulations based on using the expressions from [2] run about 200 times faster and 105 Monte Carlo trials could typically be used.

The full-scale simulations operating on 36×36 matrices were carried out using necessarily a parallel mode with eight 64-bit nodes. The code was written in FORTRAN 77 and built for Linux OS using a Portland Group FORTRAN compiler and Open MP. A 64 by 45-point grid on the unit sphere and 400×40 points in time-domain were typical. The mesh used was uniform in cosθ and [var phi]. On average, for a mesh of this size, computation time was about 15 hours. However, by using homotopy and adaptive meshes [11], one can hope for a considerable reduction in computation time. Attempts are currently underway to accomplish this.

5. Illustrative Examples

The following examples are chosen to illustrate the calculations.

  1. Fig. 3 displays 1D time domain results [5] for the dipolar interaction of d = 15 MHz (corresponding to 5.3 G), representing the distance of 15.1 Å between the two 14N nitroxides, with B1= 30 G, and B0 = 6,200 G. The Fourier transforms are also included, which shows the respective Pake doublets. The uncorrelated case was simulated using a Monte Carlo method with random angles λ1,2, η, θ, [var phi] and a Gaussian distribution in distances with FWHM of 0.75 MHz. The upper figures represent the rigorous result and the bottom figures were computed using Eqs. III.1-3 in App. III. The doublets show small peaks at 3d/2 and weak shoulders stretching up to 3d (45 MHz), caused by the pseudosecular terms in H12. The results are very close and they show that at typical conditions, e.g. with MTSL labeled proteins (where distances are typically greater than 15 Å and the spin label tether is rather flexible), virtually undistorted (by the pseudosecular term) dipolar spectra can be obtained by using a readily achievable B1 ~ 30 G.
    Figure 3
    A time-domain 1D DQC signals and their Fourier transforms for 14N nitroxides with their magnetic tensor axes orientations distributed isotropically in the molecular frame (i.e. referred to as uncorrelated case). (Top) – A computation result based ...
  2. Fig. 4 illustrates the time-domain 2D DQC signal computed rigorously for B0=6,200 G, B1 = 60 G, d = 25 MHz. In the dipolar dimension (tdip) the signal shows pronounced oscillations and in the echo dimension (techo) it produces the echo shape. Note that the echo symmetry plane is tilted due to the dipolar evolution that occurs along the echo dimension.
    Figure 4
    Time domain 2D DQC signal is shown as 3D stack plot and contour plot. The simulation were carried out rigorously for B0 = 6,200 G, B1 = 60G, d = 25 MHz and uncorrelated 14N nitroxides. The tilt of spin echo refocusing line is clearly visible. The reason ...
  3. Fig. 5 serves to illustrate the main concept of 2D FT DQC. The (filled) contour plots were produced by two-dimensional FT and are shown in the magnitude mode. Coupling between techo and tdip was removed by the shear transformation conducted in the frequency domain as fechofecho+ fdiptecho/2Δtdip) and the constant dipolar background signal (cf. Ref.3) has been removed in all 2D FT plots shown. Results for uncorrelated 14N and 15N nitroxide pairs with 2 MHz dipolar coupling (corresponding to ca. 30 Å) are shown on the left side with fixed rigid arrangement with λ1,2=(0°, 90°, 0°; 0°, 90°, 0°) shown on the right side. The B1 was infinite by using Hp [proportional, variant] Sx and the pseudosecular term in H12 was set to zero. The 2D FT spectra were summed over the range of ESR frequencies to produce a 1D dipolar spectrum on the right side of the 2D plot and over the range of dipolar frequencies to produce 1D ESR spectra at the top of each 2D plot. The Pake doublets thus correspond to a 1D FT experiment (such as shown in Fig. 3). Note that there is virtually no difference in the 1D dipolar spectra from uncorrelated and correlated cases, as one would expect for the strong pulses, which excite all orientations. However, this hidden information is developed in the 2D representation, wherein the uncorrelated case shows no variation of dipolar spectrum along the ESR dimension, whereas in the correlated case there is clearly a distinct pattern of such variations.
    Figure 5
    2D DQC (filled) magnitude contour plots obtained by 2D FT with respect to tdip and techo. Top row - 14N uncorrelated (a) and correlated (b) case. Bottom row – 15N uncorrelated (c) and correlated (d) cases. B0 = 6,200 G, d = 2 MHz. B1 was set to ...
  4. Fig. 6 serves to demonstrate the advantage of 2D spectra. The simulations were carried out rigorously for 15N nitroxides at Q-band with the following simulation parameters: B0 = 12,500 G, B1 = 60 G, d = 25 MHz has a Gaussian distribution with FWHM of 5 MHz. The angles beta were set to 90° for the correlated case (b). Even though the 1D spectrum is nearly completely smeared due to the distribution in d, the 2D plot for the correlated case is very distinct from a pattern for the uncorrelated case (a). The latter is composed of features aligned predominantly parallel to the ESR frequency axis, as one can see more clearly in Fig. 5.
    Figure 6
    2D DQC (filled) magnitude contour plots computed for 15N nitroxides using B0 = 12,500, B1 = 60 G, d = 25 MHz with a Gaussian distribution in d (FWHM = 5 MHz). Panel (a) represents an uncorrelated case, whereas in (b) the angles beta were (90°, ...
  5. Figs. 7, ,88 illustrate other cases but using two different plotting formats: the former show stack plots and the latter shows contour plots. The 2D simulations were made rigorously for B0 = 6200 G, B1 = 60 G, d = 25 MHz (12.7 Å). The top two cases represent the situation of pronounced correlations, i.e. (β1, β2) = (0°, 0°) and (90°, 90°), with the rest of the angles (α1, α2, γ1, γ2) being zero. [Note that strong effects of pseudosecular terms are clearly visible over a broad range of distances up to at least 30−35 Å.] The sensitivity to the αk and γk angles mainly depends on their particular combination with the βk, but in most cases the former angles have much less effect. For example, if for the case of (90°, 90°) one sets α2 to 90°, this will bring it close to the (0°, 90°) case (c), where in a 1D representation it is hard to see a difference from the uncorrelated case (d). In a 1D plot (especially if distances are distributed) the Pake doublets for the case (c) (0°, 90°) and (d) (uncorrelated case) are hard, if at all possible, to distinguish; but they are quite different in the 2D plots. Note the shapes of the 1D summed spectra. The case of (a) (and less so, case (b)) clearly distinguishes itself in the 1D spectrum from the uncorrelated case (d) and mainly due to the pseudosecular terms. It means that at distances of 35 Å or greater all 1D spectra will be almost identical as in the limiting case of Fig. 5. Even relatively narrow distributions in distance would likely obscure the information on orientations in a 1D plot. The matter about the smearing effect of distance distributions also applies to the cases (a, b), however in the 2D plot the overall pattern is more immune to this, since for the 2D plot the pattern is distinctly structured and spatially well-resolved. Fig. 6 clearly supports this observation.
    Figure 7
    Examples of 2D FT magnitude stack plots for three cases of orientational correlation: angles beta are (a) − (0°, 0°); (b) (0°, 90°); (c) (90°, 90°). The other four angles were set to zero. Case (d) ...
    Figure 8
    The same cases as shown as in Fig 7, but in a contour plot representation. The magnitude 2D signal was summed along both dimensions and is shown as the 1D ESR absorption spectrum (at the top) or Pake doublets (on the rhs).
  6. An extensive example compiled of 18 plots in Figs. 9a-c is intended to test the accuracy and applicability of efficient ways of signal computation based on analytical approximations for 1D DQC spectra (cf. Appendix III). The simulations were made for three orientations in angles beta, as indicated. All simulations were made for B0 = 6,200 G. Rows a1-a3 show rigorous simulations, rows b1-b3 are computed using expressions from [2] (cf. Eqs. 1-3 in App. III). In Figs. 8a-b rows a1, b1 correspond to B1 = 30 G, d = 52 MHz, which corresponds to 10 Å; rows a2 and b2 are for B1 = 60 G, and the same d = 52 MHz. Fig. 8c shows the case of B1 = 60 G and d = 25 MHz (12.8 Å). In the comparison of a1 with b1 one can clearly see some deviations, especially for (90°, 90°), but they are smaller in the remaining cases. In a2, b2 the differences are still visible but are small enough as to have little practical significance. Finally, for the last case (shown in Fig. 8c) results are virtually indistinguishable. From these results, we can set a “border region” criterion as B12D (or B13d) for which rigorous 1D simulations should be seriously considered.
    Figure 9Figure 9Figure 9
    a, The comparison of rigorous (a1) and approximate (b1) 1D computations. Three selected cases of nitroxides orientations in angles beta were used (0°, 0°), (0°, 90°), (90°, 90°). Simulations were made for ...
  7. Figure 10 shows the effect of increasing the amplitude of B1 on the maximum of the echo signal expressed by Eq. 21. It is seen from this figure, looking at the values obtained for B1 = 10 to 100, (and the value at infinity) for d =25 MHz and B0 = 6,200G, that it reaches half the value of the basic Hahn echo in an asymptotic manner as B1 approaches infinity. This is expected from the basic theory [2]. It should be noted however that in the case of large d, when pseudosecular terms are significant, there is a loss of 5−10% of the signal even at infinite B1, due to the increase of low-frequency content in the time-domain signal.
    Figure 10
    Normalized echo amplitude as a function of B1 and frequency. Plot of the rigorously computed maximum value of the echo signal, Eq. (21), as a function of B1 for d = 25 MHz, B0 = 6,200 G, uncorrelated case (circles) and d = 15 MHz, B0 = 34,000 G, orientations ...

6. Conclusions and future prospects

The main features and conclusions from the simulations presented here are as follows.

  1. The simulations for cases of short-distances (10−15 Å) are rigorously performed using the approaches presented here, where the full spin Hamiltonian is utilized during the pulse. Furthermore, the application of a very strong B1 field leads to clean Pake doublets in 1D that enables one to determine the dipolar (and exchange) couplings with the effects of correlations with the nitroxide magnetic tensors largely suppressed, in most cases. Then, in the 2D format one may examine the “fingerprint” and make distinctions amongst different orientations of the principal-axis systems of the magnetic tensors of the nitroxides when correlations are present.
  2. It was clearly demonstrated that the concept of increased correlation sensitivity in 2D FT spectra is indeed valid.
  3. We also demonstrate the increased DQC signal strengths obtained by performing experiments with stronger pulses.
  4. The criterion for using the approximate analytical approach vs. rigorous 1D simulations at conventional frequencies (up to Q band) was established.
  5. For all practical purposes rigorous DQC simulations should be utilized for the 2D domain, strong coupling cases, and the mm-wave range. 1D simulations based on analytic approaches are two to three orders of magnitude computationally more efficient and virtually linearly scalable in multiprocessor-systems. This makes it possible to apply them to more complex cases that include averaging over multiple parameters, data fitting, or to multi-spin systems.
  6. The pseudosecular terms are useful by providing more telling 1D and 2D dipolar data. The pseudosecular part of dipolar coupling is responsible for the spectral peaks with 3d/2 splitting when two spins resonate at sufficiently close frequencies. This depends differently on orientational correlations than that for the secular part, leading to a richer 2D spectrum.
  7. In the presence of distance distributions, the 1D spectrum becomes structureless, however the 2D spectrum does exhibit patterns that are distinct from those in the absence of correlation.
  8. With the approach adopted here, the electron spins are treated in the point dipole approximation ignoring spin-density delocalization. For distances less than about 10 Å, one should account for spin-density delocalization, which is significant e.g. on tyrosyl or flavin radicals, leading to a rhombic dipolar tensor.
  9. Relaxation effects have not been considered here. Phase relaxation can be introduced phenomenologically as in [1,2], but sufficiently fast spin-lattice relaxation does require treatment with full rigor in Liouville space. However, there exist simplified versions that can be used in Hilbert space, see e.g. Lee, Patyal, and Freed [12].


The authors are thankful to Ralph “Barry” Robinson for the help with building parallel code and to Dr. Richard H. Crepeau for useful advice. Access to Cornell University Theory Center is greatly acknowledged. This research was supported by NSERC, Canada and NIH/NCRR Grant P41RR016292, USA.



The Zeeman and hyperfine part of the nitroxide spin Hamiltonian, H0, in the irreducible spherical tensor operator (ISTO) representation takes the form [13]:


where L is the tensor rank (0 or 2); M = −L, ... L; k = 1,2 numbers the nitroxides; and [ell] denotes the reference frame where the tensors are defined, i.e. magnetic frame (gk), molecular (dipolar) frame (d), or laboratory frame (l). Aμk,L,M are the spin operators with μk referring to the kind of magnetic interaction, electron Zeeman (g), nuclear Zeeman (N), or hyperfine (A), and they are usually defined in the laboratory frame, Fμj,L,M is proportional to the ISTO of the magnetic interaction and is most conveniently defined in the g-frame. The transformation of Fμk,gkL,M to the laboratory frame yields the H0. In the high-field limit, where the non-secular terms (S±, S±Iz, S±I±, S±Im ) can be omitted, the ISTO form of the g-tensor reduces to:


Similarly, the relevant components of the hyperfine tensor are:


Also, the nuclear Zeeman term (when retained) is given by kFNk,l0,0ANk,l0,0 with


It is noted that the nuclear quadrupole term for 14N nitroxide is neglected. The second rank tensors Fμ,gk2,M are transformed in two steps: first, from the k-th nitroxide g-tensor axes to the dipolar frame, with its z-axis coincident with the vector r connecting magnetic dipoles, and then to the laboratory frame. The transformations from the dipolar frame to g-frame are defined by the Euler angles λk [equivalent] (αkkk) as shown in Fig. 2 and the transformation from the laboratory frame to the dipolar frame, which is defined by η [equivalent] (0,θ,[var phi]). The transformed tensors thus can be written as:


The Hamiltonian in the laboratory frame takes the form:


The coefficients Ck, Ak, Gk, and Bk are expressed as follows:




includes the transfrormations Dm,m2(λk) from the dipolar frame to the k-th magnetic frame. Since all the transformations were carried out numerically, the explicit expressions for Ck, Ak, and Bk were unnecessary. They are somewhat lengthy and the reader is referred to [1]. We just mention that the terms Ck, Ak, and Bk contain all anisotropies in the g and hf tensors as well as the Euler angles needed for their transformation from the respective principal-axes system to the laboratory frame. It is also noted that Ck, Gk, and Ak are real, whereas Bk is complex. Finally, it is noted that in carrying out the computations for B0 well up to Q band (35 GHz), the nuclear Zeeman term can be safely omitted and the H0 in real form can be obtained to a high accuracy based on the useful approximation described in [14].



In this appendix are given the details of how to calculate the final density matrix, ρf, from which the DQC echo signal can be calculated as given by Eqs. II.3-4. Basically, this consists of carrying out a series of transformations of the density matrix by a propagator, pulse or free evolution, and choosing the matrix elements of the density matrix after the application of a pulse on a coherence pathway. Here, the direct-product representation for S1 [multiply sign in circle] S2 [multiply sign in circle] I1 [multiply sign in circle] I2 to write the matrix elements of the various operators involved in the calculation will be used. To this end, the following details are required.

(i) Matrix representation and notation

For each electron, with spin S = 1/2, the matrix dimension is 2, whereas for each nucleus, with spin I, it is (2 I+1). Therefore, for the electron-nuclear spin coupled system of the nitroxide pair the size of the product space is N×N with N = 4(2I1+1)(2I2+1), which for 14N nitroxides is 36×36. The Zeeman basis with the basis vectors km1,m2;M1,M2 is used, where k's are the eigenvectors of the z-components of the electron and nuclear spin operators: Szm=mm and IzM=MM. Here m and M are the electronic and nuclear spin magnetic quantum numbers, respectively. Lower-case Roman letters will be used to describe the basis state in the product space, Greek letters will be used to describe the eigenvectors of H:Hα=ωαα. The Hamiltonian H is diagonalized by the unitary transformation VHV=E, where E is a diagonal matrix of eigenvalues of H and V is Hermitian adjoint of V. The columns of V are the eigenvectors of H:αk=kαVka. In the computations the matrices E and V are the outputs of the matrix diagonalization subroutine, such as JACOBI [15], used here.

(ii) Initial density matrix in product space

Using the expression for ρ0 as given by Eq. 7, the initial density matrix is expressed as


where 1I1 and 1I2 are 3×3 unit matrices in the respective nuclear-spin spaces for 14N nitroxides. A diagonal matrix of order 4×4 on the right-hand side of Eq. II.1 represents Sz [equivalent] S1z+S2z in the product space S1 [multiply sign in circle] S2 for the two nitroxide electron spins, that is (σz [multiply sign in circle] 12 + 11 [multiply sign in circle] σz)/2 = diag (1,0,0,−1), where 11 and 12 are 2×2 unit matrices in the spin spaces for electrons 1 and 2, respectively.

(iii) Transformation of the density matrix by a propagator

The propagators to be considered here are either a pulse propagator, due to a π/2, or a π pulse, or a free-evolution operator in the absence of a pulse. The effects of the various propagators as shown in Fig. 1 are calculated using Eq. II.3 below, with the appropriate Hamiltonians and their durations. The procedure for a given density matrix, ρ and Hamiltonian, H, acting during the time period, t, is described here, which can be specialized for the various propagators by appropriate substitutions, using appropriate times tk (t1,2= tp; t3,4= tDQ ; t5= tmtp; t6 = t5 + techo) and Eqs. 12-16, which describe the Hamiltonians used in the analysis of the 6-pulse DQC sequence. The transformed density matrix, ρ′, under the action of a propagator is expressed as:


Since eiHkmt=VkαeiωαtVαm the matrix elements of ρ′ in Eq. II.2 can be expressed as ρjk=eiHjmtρmneiHnkt=VjαeiωαtVmαρmnVnβeiωβtVkβ, where summation is conducted over the repeating indexes, or explicitly


with ωαβ [equivalent] ωαωβ. Equation II.3 can be written using a short-hand notation as ρ′ = Lρ. In this notation L is an operator, which is Q for a free evolution period or R for the action of a pulse. Coherence pathway selection implies retaining only those elements of ρ′, which belong to the pathways of interest, with the subsequent summation conducted over all pathways that contribute to the echo of interest. In computations this is accomplished by retaining only those matrix elements that correspond to the selected pathway, setting the rest to zero. This may be expressed as the application of a projection operator P, (which in reality does not need to be constructed). The final density matrix after application of N pulses and subsequent evolution periods is then calculated as


The product is computed for the full set of coherence pathways {pk} that contribute to the echo and the sum is then taken to be finally used in computing of Tr[ρfS+ ]

(iv) Coherence pathway selection

Subsequent to the action of a pulse propagator, of the matrix elements as calculated in (iii) above, all but those in the electronic product subspace of the density matrix ρ that correspond to selected coherence order p should be set to zero. The correspondence of ρik to p is compiled in Table 1 pertinent to the coherence pathways depicted in Fig. 1 illustrating the coherence pathways of the six-pulse DQC sequence. This selection of coherence pathways is achieved experimentally through phase cycling [2] or in computations is based on Table 1.

Table 1
Coherence pathways and respective matrix elements for two coupled spins.



Here we include for completeness the equation from [2] that was used in this work for making the comparison with 1D DQC signals produced in rigorous computations. The echo amplitude, V, is a function of tdip=2tp-tm and is given by


The time variables are defined in accordance with Fig.1 and notations in the text. F(t) has the form:


Here, A = d(1−3cos2θ) and b = −A/2, represent the secular and pseudosecular parts of the dipolar coupling; R2 = Δω2 + b2, where Δω = ω1ω2 is the difference between the Larmor frequencies ω1 and ω2 of the nitroxide's electron spins in the frame of reference rotating with the frequency ωrf of the excitation pulses, which was set to coincide with the center of the nitroxide ESR spectrum. Also, q = b / R and p2 = 1− q2. The amplitude factors were taken in the simplest possible form as


where ω1rf = γeB1 for π pulses, which are taken equal here, but do not have to be.

Since ωk=ωk (λk,η), the powder averaging is conducted essentially in the same way as in rigorous computations with ωk determined for each set of (M1, M2). This was accomplished using an approximation [14].


1. Saxena S, Freed JH. J. Chem. Phys. 1997;107:1317–1340.
2. Borbat P, Freed JH. In: Biological Magnetic Resonance. Berliner LJ, Eaton SS, Eaton GR, editors. Vol. 19. Kluwer Academic/Plenum Publications; New York: 2000. pp. 383–459.
3. Borbat PP, Freed JH. Multiple-quantum ESR and distance measurements. Chem. Phys. Lett. 1999;313:145–154.
4. Borbat PP, Freed JH. Measuring distances by pulsed dipolar ESR spectroscopy: Spin-labeled histidine kinases. Methods in Enzymology. 2007;423:52–116. [PubMed]
5. Borbat PP, Freed JH. “Pros and Cons of Pulse Dipolar ESR: DQC & DEER” EPR Newsletter 17. 2007;(2−3):21–29.
6. Lee S, Budil DE, Freed JH. Theory of two-dimensional Fourier transform ESR for ordered and viscous fluids. J. Chem. Phys. 1994;99:7098–7107.
7. Suzuki M. Decomposition formulas of exponential operators and lie exponentials with some applications to quantum-mechanics and statistical physics. J. Math. Phys. 1985;26:601–612.
8. Trotter HF. On the product of semi-groups of operators. Am. Math. Soc. 1959;10:545–551. Proc.
9. Ernst RR, Bodenhausen G, Wokaun A. Principles of nuclear magnetic resonance in one and two dimensions. Clarendon Press; Oxford: 1987.
10. Gemperle C, Aebli G, Schweiger A, Ernst RR. Phase cycling in pulse EPR. J. Mag. Res. 1990;88:241–256.
11. Gates KE, Griffin M, Hanson GR, Burrage K. Computer simulation of magnetic resonance spectra employing homotopy. J. Magn. Reson. 1998;135:103–112. [PubMed]
12. Lee SY, Budil DE, Freed JH. Theory of 2-Dimensional Fourier-Transform Electron-Spin-Resonance for Ordered and Viscous Fluids. J. Chem. Phys. 1994;101:5529–5558.
13. a. Schneider DJ, Freed JH. In: Spin Labeling: Theory and Application. Berliner LJ, editor. Academic; New York: 1976. Chap. 2. b. Schneider DJ, Freed JH. In: Advances in Chemical Physics. Hirschfelder JO, Wyatt RE, Coalson RD, editors. LXIII. Wiley; New York: 1989. c. Misra SK. Simulation of slow-motion CW EPR spectrum using stochastic Liouville equation for an electron spin coupled to two nuclei with arbitrary spins: Matrix elements of the Liouville superoperator. J. Magn. Reson. 2007;189:59. [PubMed]
14. Libertini LJ, Griffith OH. Orientation dependence of the electron spin resonance spectrum of di-t-butyl nitroxide. J. Chem. Phys. 1970;53:1359–1367.
15. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in Fortran. Second Edition. Cambridge University Press; 1992. pp. 456–462. A better version of the JACOBI subroutine than that given in the book can be found on the Netlib website or else is available from the authors.