Analytical expressions for the time evolution of spin systems affected by two or more interactions

Analytical expressions for the description of the time evolution of spin systems beyond product–operator formalism (POF) can be obtained if a low-dimensional subspace of the Liouville space has been found in which the time evolution of the spin system takes place completely. This can be achieved using a procedure that consists of repeated application of the commutator of the Hamiltonian with the density operator. This iteration continues as long as the result of such a commutator operation contains a term that is linearly independent of all the operators appearing in the previous commutator operations. The coefficients of the resulting system of commutator relations can be immediately inserted into the generic propagation formulae given in this article if the system contains two, three, or four equations. In cases where the validity conditions of any of these propagation formulae are not met, the coefficients are used as intermediate steps to obtain both the Liouvillian and propagator matrices of the system. Several application examples are given where an analytical equation can be obtained for the description of the time evolution of small spin systems under the influence of two or more interactions. This procedure for finding the Liouvillian matrix is not limited to time-independent interactions. Some examples illustrate the treatment of time-dependent problems using this method.

1 Introduction

Considerable progress has been made in predicting the time evolution of spin systems. Numerical calculations or simulations are possible in large systems of coupled spins (Kuprov et al., 2007), and the evolution during an arbitrarily long sequence of pulses can be simulated with the help of software (Veshtort and Griffin, 2006; Bak et al., 2000; Hogben et al., 2011), including the effect of thermal motion.

However, there may be situations where the analytical representation of the time evolution is advantageous, providing more physical intuition than a numerical procedure, which can seem like a black box. The desire or need to take a look inside this black box could be a motivation to deal with analytical contexts. In many cases, product–operator formalism (POF) (Slichter, 1987; Sorensen et al., 1983; Packer and Wright, 1983; Wang and Slichter, 1986) has been used for this purpose. This is a rather intuitive scheme where the states of the spin system are represented by spin operators. Such a tool can be very useful for doing short calculations without access to fundamental quantum mechanics or without a simulation program. For demonstrations in discussions and lectures, e.g., on the effects of pulse sequences like INEPT or HSQC (Van De Ven and Hilbers, 1983), an illustrative explanation can be given here.

As an example, consider the POF description of the propagation of transversal magnetization of spin I = $\mathrm{1}/\mathrm{2}$ which is coupled by a scalar interaction to another spin S = $\mathrm{1}/\mathrm{2}$ (coupling constant J):

$$\begin{array}{}\text{(1)}& {\widehat{I}}_x\stackrel{\mathrm{2}\mathit{\pi }J{\widehat{I}}_z{\widehat{S}}_z\cdot t}{⟶}{\widehat{I}}_x\cos \mathit{\pi }Jt+\mathrm{2}{\widehat{I}}_y{\widehat{S}}_z\sin \mathit{\pi }Jt.\end{array}$$

This means that the spin system oscillates between two states characterized by the operators ${\widehat{I}}_z$ and $\mathrm{2}{\widehat{I}}_z{\widehat{S}}_z$. Thus, the evolution takes place in a 2-D subspace of the total operator space. In contrast, to describe a system with two spins $\mathrm{1}/\mathrm{2}$, we need a 4×4 density matrix. This means that the Liouville–von Neumann equation is a system of 16 scalar differential equations for each of the matrix elements. However, obviously it is possible in this example to reduce the 16-D problem to the 2-D one described by Eq. (1).

A question arises: are there other situations where dimensionality reduction is possible? Candidates for such situations are, for example, the evolution of a spin system under dipoledipole interaction and simultaneous rf irradiation, or cross-polarization with respect to the finite RF power.

$$\begin{array}{c}{\displaystyle }{\widehat{I}}_z\stackrel{\text{dipole–dipole interaction}+\text{rf irradiation}}{⟶}\mathrm{?}\\ {\displaystyle }{\widehat{I}}_x\stackrel{\text{cross-polarization under a finite rf power}}{⟶}\mathrm{?}\end{array}$$

Are there others? For numerical calculations, any reduction in dimensionality leads to a reduction in computation time, but to enable analytical calculations, which is the goal here, this may be essential. The aims of this paper are (i) to introduce a procedure for reducing the dimension of the problem, (ii) to show examples of the application of this procedure to obtain analytical equations, and (iii) to show the resulting generic propagation rules as templates for cases where the problem could be reduced to a 2-D, 3-D, or 4-D problem. The latter can be seen as an extension of the product–operator formalism to somewhat more complex situations. (i) and (ii) can be useful for simplifying some calculations through dimension reduction without any approximation. Even for numerical calculations, this can be useful if it helps us to work on low-dimensional systems. As an example, a system of two spins I and S is mentioned here, which is subject to dipolar interaction but which is decoupled by rf irradiation, and at the same time magic-angle sample spinning (MAS) modulating the dipolar oscillation takes place. Its time evolution can be described by a system of three differential equations using the method presented here. The application of the Shirley–Floquet method will be greatly simplified here. On the other hand, the application of the Liouville–von Neumann equation in the 4-D wave-function space leads to a system of 16 differential equations, even if some of its coefficients can be 0.

In Sect. 2, the mathematical background is investigated, which makes it possible to reduce the dimension to the value 2 in cases where the POF is applicable. The results of this calculation are applied to more complex structures in Sect. 3. Finally, in Sect. 4, some template formulae are given, together with examples of cases where a reduction to 3-D, 4-D, 5-D, and 6-D problems is possible.

As usual, in this paper operators are denoted by a hat ($\widehat{A}$) and superoperators by a double hat ($\widehat{\widehat{B}}$), while the vector or matrix associated with an operator is denoted by the same symbol but in italics and in bold roman, respectively, and without a hat (A and B). Scalar variables are written in italics.

2 Dimension reduction through POF

2.1 Condition for validity: commutator relations

As shown in the references cited above, the generic scheme of the POF is as follows. The time evolution can be predicted by the propagation rules

$$\begin{array}{}\text{(2)}& \begin{array}{l}\widehat{A}\stackrel{\widehat{H}t}{⟶}\widehat{A}\cos \mathit{\lambda }t+\widehat{B}\sin \mathit{\lambda }t\text{and}\\ \widehat{B}\stackrel{\widehat{H}t}{⟶}\widehat{B}\cos \mathit{\lambda }t-\widehat{A}\sin \mathit{\lambda }t\end{array}\end{array}$$

if and only if

$$\begin{array}{}\text{(3)}& \begin{array}{l}\left[\widehat{H},\widehat{A}\right]=i\mathit{\lambda }\widehat{B}\text{and}\\ \left[\widehat{H},\widehat{B}\right]=-i\mathit{\lambda }\widehat{A},\end{array}\end{array}$$

where $\left[\widehat{A},\widehat{B}\right]\equiv \widehat{A}\widehat{B}-\widehat{B}\widehat{A}$ denotes the commutator between the operators $\widehat{A}$ and $\widehat{B}$. That is, Eq. (2) describes the motion of the density operator in a 2-D subspace of the total Liouville space of the current spin system, although the Liouville space has a much larger dimension.

However, Eq. (3) is often not satisfied when more than one interaction must be considered. This applies, for example, to an rf irradiation with strength ω₁ in an ensemble of spins I = $\mathrm{1}/\mathrm{2}$ which are coupled to another spin ensemble S = $\mathrm{1}/\mathrm{2}$ (coupling frequency D_IS). This situation can be described by the Hamiltonian $\widehat{H}=-\mathrm{2}{D}_{\mathrm{IS}}{\widehat{I}}_z{\widehat{S}}_z-{\mathit{\omega }}_{\mathrm{1}}{\widehat{I}}_x$ and an initial state ${\mathit{\rho }}_{\mathrm{0}}={\widehat{I}}_z$. The double calculation of the commutator of the Hamiltonian with ${\widehat{I}}_z$ results not only in ${\widehat{I}}_z$, but also an additional term. In this case, after the third application of this commutator operation, the result will contain only the operators already used:

$$\begin{array}{}\text{(4)}& \begin{array}{l}\left[\widehat{H},{\widehat{I}}_z\right]=i{\mathit{\omega }}_{\mathrm{1}}{\widehat{I}}_y,\\ \left[\widehat{H},{\widehat{I}}_y\right]=-i{D}_{\mathrm{IS}}\cdot \mathrm{2}{\widehat{I}}_x{\widehat{S}}_z-i{\mathit{\omega }}_{\mathrm{1}}{\widehat{I}}_z,\\ \left[\widehat{H},\mathrm{2}{\widehat{I}}_x{\widehat{S}}_z\right]=i{D}_{\mathrm{IS}}\cdot {\widehat{I}}_y.\end{array}\end{array}$$

The connection between the commutator Eq. (3) and the propagation Eq. (2) becomes much clearer if we reformulate both sets of equations in matrix form,

$$\begin{array}{}\text{(5)}& \left(\begin{array}{c}\widehat{A}\\ \widehat{B}\end{array}\right)\stackrel{\widehat{H}\cdot t}{⟶}\left(\begin{array}{cc}\cos \mathit{\lambda }t& \sin \mathit{\lambda }t\\ -\sin \mathit{\lambda }t& \cos \mathit{\lambda }t\end{array}\right)\left(\begin{array}{c}\widehat{A}\\ \widehat{B}\end{array}\right),\end{array}$$

if and only if

$$\begin{array}{}\text{(6)}& \left(\begin{array}{c}\left[\widehat{H},\widehat{A}\right]\\ \left[\widehat{H},\widehat{B}\right]\end{array}\right)=\left(\begin{array}{cc}\mathrm{0}& i\mathit{\lambda }\\ -i\mathit{\lambda }& \mathrm{0}\end{array}\right)\left(\begin{array}{c}\widehat{A}\\ \widehat{B}\end{array}\right).\end{array}$$

For this analysis, it is important to note that the 2×2 matrix in Eq. (5) is the exponential of the 2×2 matrix in Eq. (6) multiplied by −it (see Sect. S2 in the Supplement):

$$\begin{array}{}\text{(7)}& \left(\begin{array}{cc}\cos \mathit{\lambda }t& \sin \mathit{\lambda }t\\ -\sin \mathit{\lambda }t& \cos \mathit{\lambda }t\end{array}\right)=\exp \left(\begin{array}{cc}\mathrm{0}& \mathit{\lambda }t\\ -\mathit{\lambda }t& \mathrm{0}\end{array}\right).\end{array}$$

Obviously, for this case, it is possible to take a particular parameter (λ) of a system of commutator equations and insert it into the template Eq. (2), which is formulated as a set of two propagation formulae. Then the following questions arise:

• 1.

  Is it possible to apply this procedure to more complex cases, such as that of Eq. (4)?

• 2.

  To what dimension can we reduce a given problem?

• 3.

  Are there any template propagation formulae for the case where a reduction to a 3-D or 4-D case is possible?

As a final remark in this subsection, we note that we have performed these calculations in the operator space (Liouville space). This seems advantageous because it is the natural way of treating this topic. We describe states in terms of operators rather than wave functions (see above). Some properties of the Liouville space that are important for this work are listed in Sect. 2.2.

2.2 Properties of the Liouville space that are important in this article

In some papers, the space of wave functions is denoted as the Hilbert space in order to have a contrast with the Liouville space. However, this does not seem appropriate since the space of linear operators also satisfies the conditions of being a Hilbert space (Jordan, 2005).

• Definition of the scalar product of two operators $\widehat{A}$ and $\widehat{B}$: $\left(\widehat{A},\widehat{B}\right):=\text{Tr}\left({\mathbf{A}}^{†}\cdot \mathbf{B}\right)$, where A and B are the matrix representations of $\widehat{A}$ and $\widehat{B}$ in the wave-function space and the superscript † means the Hermitian conjugate of the corresponding matrix.

• Two operators $\widehat{A}$ and $\widehat{B}$ are said to be orthogonal if their scalar product is 0: $\widehat{A}⟂\widehat{B}\iff \left(\widehat{A},\widehat{B}\right)=\mathrm{0}$.

• The Euclidean norm of an operator is defined as the square root of the scalar product of the operator with itself: $\parallel \widehat{A}\parallel =\sqrt{\left(\widehat{A},\widehat{A}\right)}=\sqrt{\text{Tr}\left({\mathbf{A}}^{†}\cdot \mathbf{A}\right)}$. For the sake of brevity, the term “norm” is used throughout this paper to refer to the Euclidean norm. The rules for calculating the norm values for some of the operators used in this article are given in Appendix A. The norm of an operator is invariant with respect to a unitary transformation, but it depends on the relevant space, which is different for different numbers of spins. Therefore, the table in Appendix A contains different norm equations for the same operator but in different spaces.

• Each operator of this space can be expanded into a series of basis operators. In particular, for the density operator ρ, we have

  $$\begin{array}{}\text{(8)}& \widehat{\mathit{\rho }}={\mathit{\rho }}_{\mathrm{1}}{\widehat{u}}_{\mathrm{1}}+\mathrm{\dots }+{\mathit{\rho }}_d{\widehat{u}}_d={\mathit{u}}^{\mathrm{T}}\cdot \mathit{\rho },\end{array}$$

  where d is the dimension of the particular Liouville space, ρ is a column matrix which contains the expansion coefficients ρ_i (i∈{1…d}) of the density operator, u is a column matrix whose elements are the basis operators, and A^T is the transpose of the matrix A. If all the basis operators are pairwise orthogonal, ρ_i can be calculated as follows:

  $$\begin{array}{}\text{(9)}& {\mathit{\rho }}_i=\left(\widehat{\mathit{\rho }},{\widehat{u}}_i\right)/\parallel {\widehat{u}}_i\parallel \left(i\in \mathit{\left\{}\mathrm{1}\mathrm{\dots }d\mathit{\right\}}\right).\end{array}$$

• Mappings between operators are described by superoperators:

  • i.

    Liouville superoperator (or simply “Liouvillian”) for forming the commutator with the Hamiltonian:
    $\widehat{\widehat{L}}$: $\widehat{A}\mapsto \left[\widehat{H},\widehat{A}\right]$, also written as $\widehat{\widehat{L}}\widehat{A}=\left[\widehat{H},\widehat{A}\right]$,
    and

  • ii.

    propagation superoperator (or simply “superpropagator”) for describing the time evolution:
    $\widehat{\widehat{U}}$: $\widehat{\mathit{\rho }}\left(\mathrm{0}\right)\mapsto \widehat{\mathit{\rho }}\left(t\right)$, also written as

    $$\begin{array}{}\text{(10)}& \widehat{\mathit{\rho }}\left(t\right)=\widehat{\widehat{U}}\widehat{\mathit{\rho }}\left(\mathrm{0}\right).\end{array}$$

• The time evolution of the density operator is governed by the Liouville–von Neumann equation. It is formulated in the Liouville space as

  $$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\widehat{\mathit{\rho }}=-i\widehat{\widehat{L}}\widehat{\mathit{\rho }}.\end{array}$$

  If $\widehat{\widehat{L}}$ does not depend on time, the formal solution is

  $$\begin{array}{}\text{(12)}& \widehat{\mathit{\rho }}\left(t\right)=\widehat{\widehat{U}}\left(t\right)\widehat{\mathit{\rho }}\left(\mathrm{0}\right)\text{with}\widehat{\widehat{U}}\left(t\right)=\exp \left(-i\widehat{\widehat{L}}t\right).\end{array}$$

  Here we find a similarity to the matrix formulation of the POF (Eqs. 5–7). In fact, the coefficient matrix of the system of commutator equations is the transposition of the Liouvillian matrix, and the coefficient matrix of the POF (Eq. 5) is the transposition of the superpropagator. This is proven in the Supplement (Sect. S1.1).

• The norm of the density matrix is time-invariant. This can be proven by multiplying the Liouville–von Neumann Eq. (11) by $\widehat{\mathit{\rho }}$:

  $$\begin{array}{}\text{(13)}& \left(\widehat{\mathit{\rho }}.{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\widehat{\mathit{\rho }}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(\widehat{\mathit{\rho }}.\widehat{\mathit{\rho }}\right)=-i\left(\widehat{\mathit{\rho }}.\widehat{\widehat{L}}\widehat{\mathit{\rho }}\right)=\mathrm{0}.\end{array}$$

  The scalar product of a Hermitian operator ($\widehat{\mathit{\rho }}$) and its commutator with another Hermitian operator ($\widehat{H}$) is 0 (Sect. S1.2).

The Liouville space formulation objective of this article is to find a subspace of the total Liouville space that contains all possible states occurring in the current problem and that has the smallest possible dimension. An analogous procedure is not possible in the wave-function space, where the time evolution has to be calculated using $\widehat{U}\left(t\right)\widehat{\mathit{\rho }}\left(\mathrm{0}\right){\widehat{U}}^{-\mathrm{1}}\left(t\right)$. In the following section, commutator equations representing the action of the Liouvillian are established. Their coefficients are needed for further calculation of the propagation formulae.

2.3 Two forms and two bases for symbolic description of the time evolution

For practical calculations, Eq. (10) is usually transformed into a matrix equation. In the literature this is usually done in two different forms, both of which can be found, e.g., in Ernst et al. (1987) (Sect. 2.1.4).

Form 1:

$$\begin{array}{}\text{(14)}& \mathit{\rho }\left(t\right)=\mathbf{U}\left(t\right)\cdot \mathit{\rho }\left(\mathrm{0}\right).\end{array}$$

The time evolution is described here in the space of N-dimensional column matrices: ρ(t) is the density matrix according to Eq. (8), i.e., a column containing the coefficients of the expansion of the density operator $\widehat{\mathit{\rho }}\left(t\right)$ into basis operators. That is, we follow the propagation into the space ℂ^N, i.e., the space of N-row column matrices with the basis $\mathit{\left\{}\right(\mathrm{1},\mathrm{0},\mathrm{\dots },\mathrm{0}{)}^{\mathrm{T}},(\mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{0}{)}^{\mathrm{T}},\mathrm{\dots },(\mathrm{0},\mathrm{0},\mathrm{\dots },\mathrm{1}{)}^{\mathrm{T}}\mathit{\}}$. The superpropagator $\widehat{\widehat{U}}$ is an isomorphism in this space and can be represented by the N×N matrix U.

Form 2 (denoted here as the “propagation formula”):

$$\begin{array}{}\text{(15)}& \begin{array}{l}{\widehat{A}}_{\mathrm{1}}\stackrel{\widehat{H}t}{⟶}{V}_{\mathrm{11}}{\widehat{A}}_{\mathrm{1}}+V_{\mathrm{12}}{\widehat{A}}_{\mathrm{2}}+\mathrm{\dots }+V_{\mathrm{1}N}{\widehat{A}}_N,\\ {\widehat{A}}_{\mathrm{2}}\stackrel{\widehat{H}t}{⟶}{V}_{\mathrm{21}}{\widehat{A}}_{\mathrm{1}}+V_{\mathrm{22}}{\widehat{A}}_{\mathrm{2}}+\mathrm{\dots }+V_{\mathrm{2}N}{\widehat{A}}_N,\\ \mathrm{\dots }\\ {\widehat{A}}_N\stackrel{\widehat{H}t}{⟶}{V}_{N\mathrm{1}}{\widehat{A}}_{\mathrm{1}}+V_{N\mathrm{2}}{\widehat{A}}_{\mathrm{2}}+\mathrm{\dots }+V_{NN}{\widehat{A}}_N.\end{array}\end{array}$$

This means that the density operator, which is ${\widehat{A}}_{\mathrm{1}}$ at t = 0, evolves under the influence of the Hamiltonian $\widehat{H}$ during time t into a linear combination of independent operators ${\widehat{A}}_{\mathrm{1}},{\widehat{A}}_{\mathrm{2}},\mathrm{\dots },{\widehat{A}}_N$. Unlike form 1, here we work in the operator space with the basis $\mathit{\{}{\widehat{A}}_{\mathrm{1}},{\widehat{A}}_{\mathrm{2}},\mathrm{\dots },{\widehat{A}}_N\mathit{\}}$.

Connection between both forms. U=V^T if the basis column matrices (form 1) are assigned to the corresponding basis operators of the second basis. For proof, see Sect. S1.1.

This means that we can rewrite Eq. (15) as

$$\begin{array}{}\text{(16)}& \begin{array}{l}{\widehat{A}}_{\mathrm{1}}\stackrel{\widehat{H}t}{⟶}{U}_{\mathrm{11}}{\widehat{A}}_{\mathrm{1}}+U_{\mathrm{21}}{\widehat{A}}_{\mathrm{2}}+\mathrm{\dots }+U_{N\mathrm{1}}{\widehat{A}}_N,\\ {\widehat{A}}_{\mathrm{2}}\stackrel{\widehat{H}t}{⟶}{U}_{\mathrm{12}}{\widehat{A}}_{\mathrm{1}}+U_{\mathrm{22}}{\widehat{A}}_{\mathrm{2}}+\mathrm{\dots }+U_{N\mathrm{2}}{\widehat{A}}_N,\\ \mathrm{\dots }\\ {\widehat{A}}_N\stackrel{\widehat{H}t}{⟶}{U}_{\mathrm{1}N}{\widehat{A}}_{\mathrm{1}}+U_{\mathrm{2}N}{\widehat{A}}_{\mathrm{2}}+\mathrm{\dots }+U_{NN}{\widehat{A}}_N.\end{array}\end{array}$$

The set of all U values is the dual space of the set of all V. Similarly, the coefficient matrix of the commutator equations is the transposed Liouvillian matrix (see Sect. S1.1).

3 Procedure for finding propagation formulae

3.1 Requirements for a suitable subspace

To be sure that a given subspace 𝕊 contains the whole evolution of a spin system, we have to check that the action of the propagator $\widehat{\widehat{U}}$ on each operator $\widehat{A}$ of this subspace also results in an element of the subspace:

$$\begin{array}{}\text{(17)}& \forall \widehat{A}\in \mathbb{S}:\widehat{\widehat{U}}\widehat{A}\in \mathbb{S}.\end{array}$$

A subspace with this property is said to be propagator-invariant (Jordan, 2005). This is equivalent to the Liouvillian invariance of the subspace, i.e.,

$$\begin{array}{}\text{(18)}& \forall \widehat{A}\in \mathbb{S}:\widehat{\widehat{L}}\widehat{A}\in \mathbb{S},\end{array}$$

because the propagator is the sum of repeated $\widehat{\widehat{L}}$ actions:

$$\begin{array}{}\text{(19)}& \widehat{\widehat{U}}\widehat{A}=\exp \left(-i\widehat{\widehat{L}}t\right)\widehat{A}=\sum _{n=\mathrm{0}}^{\mathrm{\infty }}{\displaystyle \frac{(-it{)}^n}{n\mathrm{!}}}{\widehat{\widehat{L}}}^n\widehat{A}\in \mathbb{S}.\end{array}$$

A sufficient condition for a subspace spanned by ${\widehat{u}}_{\mathrm{1}}\mathrm{\dots }{\widehat{u}}_N$ (Eq. 8) to be Liouvillian-invariant is that the action of the Liouvillian of any of the ${\widehat{u}}_i$ values also results in an element of that subspace, i.e., $\widehat{\widehat{L}}{\widehat{u}}_i\in \mathbb{S}\forall i\in \mathit{\left\{}\mathrm{1}\mathrm{\dots }N\mathit{\right\}}$. A Liouvillian-invariant subspace which is of interest here should at least contain the initial density operator. Furthermore, all multiple actions of the Liouvillian on the density operator must result in elements of that subspace: ${\widehat{\widehat{L}}}^n\widehat{\mathit{\rho }}\left(\mathrm{0}\right)\in \mathbb{S}\forall n\in \mathbb{N}$.

In principle we can construct such a subspace as the set of all linear combinations of $\widehat{\mathit{\rho }}\left(\mathrm{0}\right)$, $\widehat{\widehat{L}}\widehat{\mathit{\rho }}\left(\mathrm{0}\right)$, ${\widehat{\widehat{L}}}^{\mathrm{2}}\widehat{\mathit{\rho }}\left(\mathrm{0}\right),\mathrm{\dots },{\widehat{\widehat{L}}}^N\widehat{\mathit{\rho }}\left(\mathrm{0}\right)$, where N is the largest number for which this operator set is linearly independent. This means that ${\widehat{\widehat{L}}}^{N+\mathrm{1}}\widehat{\mathit{\rho }}\left(\mathrm{0}\right)$ can be represented as a linear combination of all lower powers. Then, all further applications of $\widehat{\widehat{L}}$ lead to operators which are also linearly dependent.

This procedure is reminiscent of the formation of Krylov subspaces in matrix spaces (Watkins, 2007). For further considerations in this article, in particular the calculation of matrices, it is convenient to have basis operators that are pairwise orthogonal. This can be achieved by combining the Krylov-like procedure with the Gram–Schmidt orthogonalization, which is analogous to the Arnoldi procedure for matrix spaces (Watkins, 2007). The details of the procedure are presented in the next subsection.

3.2 First step: creating a closed system of commutator equations

According to the result of the previous subsection, the application of the Liouvillian to the initial density operator is repeated as long as the result contains a component that is linearly independent of all previous results. The action of the Liouvillian consists in forming the commutator of the Hamiltonian with the considered operator. Consequently, the search for this Liouvillian-invariant minimum subspace is performed by repeatedly calculating the commutator of the Hamiltonian using the operator representing the initial state of the spin system, as shown in detail below.

Let us assume that the system under consideration is characterized by the Hamiltonian $\widehat{H}$ and the initial state density operator by ${\widehat{A}}_{\mathrm{1}}$.

Evaluation of the first commutator: $\left[\widehat{H},{\widehat{A}}_{\mathrm{1}}\right]$

If the result is 0, the state of the spin system is constant in time; see the Supplement (examples 0D-1 and 0D-2). In the case where this commutator does not vanish, it can be decomposed into a term which is proportional to ${\widehat{A}}_{\mathrm{1}}$ and another term ${\widehat{A}}_{\mathrm{12}}$ which is orthogonal to the first one:

$$\begin{array}{}\text{(20)}& \left[\widehat{H},{\widehat{A}}_{\mathrm{1}}\right]={\mathit{\lambda }}_{\mathrm{11}}{\widehat{A}}_{\mathrm{1}}+{\widehat{A}}_{\mathrm{12}}\text{with}{\widehat{A}}_{\mathrm{12}}⟂{\widehat{A}}_{\mathrm{1}}.\end{array}$$

The commutator of two Hermitian operators is always orthogonal to both, as proven in the Supplement. This means that, for a Hermitian ${\widehat{A}}_{\mathrm{1}}$, $\left[\widehat{H},{\widehat{A}}_{\mathrm{1}}\right]$ cannot have a component that contains ${\widehat{A}}_{\mathrm{1}}$ itself. This is different for non-Hermitian operators (see Sect. 4.1).

We replace ${\widehat{A}}_{\mathrm{12}}\to {\mathit{\lambda }}_{\mathrm{12}}{\widehat{A}}_{\mathrm{2}}$, choosing the scalar λ₁₂ so that ${\widehat{A}}_{\mathrm{2}}$ has the same norm as ${\widehat{A}}_{\mathrm{1}}$. The corresponding table in Appendix A can be used to determine the norms of the operators.

Evaluation of the second commutator: $\left[\widehat{H},{\widehat{A}}_{\mathrm{2}}\right]$

This will be decomposed as

$$\begin{array}{}\text{(21)}& \left[\widehat{H},{\widehat{A}}_{\mathrm{2}}\right]={\mathit{\lambda }}_{\mathrm{21}}{\widehat{A}}_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{22}}{\widehat{A}}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{23}}{\widehat{A}}_{\mathrm{3}},\end{array}$$

with the condition that ${\widehat{A}}_{\mathrm{3}}$ is orthogonal to both ${\widehat{A}}_{\mathrm{1}}$ and ${\widehat{A}}_{\mathrm{2}}$. The coefficient λ₂₃ is chosen so that ${\widehat{A}}_{\mathrm{3}}$ has the same norm as ${\widehat{A}}_{\mathrm{1}}$ and ${\widehat{A}}_{\mathrm{2}}$.

If ${\widehat{A}}_{\mathrm{3}}$ = 0, i.e., $\left[\widehat{H},{\widehat{A}}_{\mathrm{2}}\right]$ is a linear combination of ${\widehat{A}}_{\mathrm{1}}$ and ${\widehat{A}}_{\mathrm{2}}$, the procedure is finished. Then we have a system of commutator equations like Eq. (3), i.e., the usual POF.

Evaluation of the nth commutator: $\left[\widehat{H},{\widehat{A}}_n\right]$

Similarly, it will be expanded into a series of operators known from the previous commutator evaluations and a remainder:

$$\begin{array}{}\text{(22)}& \left[\widehat{H},{\widehat{A}}_n\right]={\mathit{\lambda }}_{n\mathrm{1}}{\widehat{A}}_{\mathrm{1}}+\mathrm{\dots }+{\mathit{\lambda }}_{n,n}{\widehat{A}}_n+{\mathit{\lambda }}_{n,n+\mathrm{1}}{\widehat{A}}_{n+\mathrm{1}},\end{array}$$

with the condition that ${\widehat{A}}_{n+\mathrm{1}}$ is orthogonal to all of the ${\widehat{A}}_k$ values with k∈1…n. Again, the coefficient ${\mathit{\lambda }}_{n,n+\mathrm{1}}$ is chosen so that ${\widehat{A}}_{n+\mathrm{1}}$ has the same norm as the other operators of this set.

End of the procedure. If for a certain n = N the commutator is a linear combination of the previously determined ${\widehat{A}}_k$ without any remainder, the iteration is finished. The set of pairwise orthogonal operators $\mathit{\left\{}{\widehat{A}}_{\mathrm{1}}\mathrm{\dots }{\widehat{A}}_N\mathit{\right\}}$ spans a Liouvillian-invariant subspace of the entire Liouville space. Its dimension is N.

Possible modification of the procedure. The remainder of any commutator evaluation can be written as the sum of two or more operators, all of which must be orthogonal to the other operators. In some cases, this can simplify the coefficient matrix.

3.3 Second step: Liouvillian matrix

The result of the whole procedure is a system of equations like the following:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}\left(\begin{array}{c}\left[\widehat{H},{\widehat{A}}_{\mathrm{1}}\right]\\ \left[\widehat{H},{\widehat{A}}_{\mathrm{2}}\right]\\ \left[\widehat{H},{\widehat{A}}_{\mathrm{3}}\right]\\ \mathrm{\dots }\\ \left[\widehat{H},{\widehat{A}}_{N-\mathrm{2}}\right]\\ \left[\widehat{H},{\widehat{A}}_{N-\mathrm{1}}\right]\\ \left[\widehat{H},{\widehat{A}}_N\right]\end{array}\right)=& \left(\begin{array}{ccccccc}{\mathit{\lambda }}_{\mathrm{11}}& {\mathit{\lambda }}_{\mathrm{12}}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\mathit{\lambda }}_{\mathrm{21}}& {\mathit{\lambda }}_{\mathrm{22}}& {\mathit{\lambda }}_{\mathrm{23}}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathit{\lambda }}_{\mathrm{32}}& {\mathit{\lambda }}_{\mathrm{33}}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& {\mathit{\lambda }}_{N-\mathrm{2},N-\mathrm{2}}& {\mathit{\lambda }}_{N-\mathrm{2},N-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& {\mathit{\lambda }}_{N-\mathrm{1},N-\mathrm{2}}& {\mathit{\lambda }}_{N-\mathrm{1},N-\mathrm{1}}& {\mathit{\lambda }}_{N-\mathrm{1},N}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& {\mathit{\lambda }}_{N,N-\mathrm{1}}& {\mathit{\lambda }}_{N,N}\end{array}\right)\\ & \left(\begin{array}{c}{\widehat{A}}_{\mathrm{1}}\\ {\widehat{A}}_{\mathrm{2}}\\ {\widehat{A}}_{\mathrm{3}}\\ \mathrm{\dots }\\ {\widehat{A}}_{N-\mathrm{2}}\\ {\widehat{A}}_{N-\mathrm{1}}\\ {\widehat{A}}_N\end{array}\right).\end{array}\end{array}$$

As shown in Sect. S1.1, the coefficient matrix is the transposition of the Liouvillian matrix. This system of equations implies that the action of the Liouvillian on any ${\widehat{A}}_i(i\in \mathit{\{}\mathrm{1}\mathrm{\dots }N\mathit{\left\}}\right)$ leads to a linear combination of ${\widehat{A}}_i$. In other words, the subspace spanned by the operators ${\widehat{A}}_{\mathrm{1}}\mathrm{\dots }{\widehat{A}}_n$ is both Liouvillian-invariant and propagator-invariant. The density operator, once located in this subspace, will not leave it as long as the interaction does not vary. This explains why the POF can be applied successfully as a 2-D problem, even if the complete Liouville space has a much higher dimension.

The zeros in the upper triangle of the coefficient matrix result from the above: if the remainder of the nth commutator is identified with only one new operator, then ${\mathit{\lambda }}_{n,n+\mathrm{1}}$ and ${\widehat{A}}_{n+\mathrm{1}}$ are determined. ${\mathit{\lambda }}_{n,n+\mathrm{2}},\mathrm{\dots }$ and ${\widehat{A}}_{n+\mathrm{2}},\mathrm{\dots }$ are still unknown in this step; the matrix elements to the right of them remain zero. In the case of the modified procedure mentioned above, the matrix structure in Eq. (23) changes. Then, ${\mathit{\lambda }}_{n,n+\mathrm{2}}\mathrm{\dots }$ can also be nonzero. Due to the hermiticity of the Liouvillian, i.e., $L_{ij}=L_{ji}^{*}$, the lower triangle matrix has the same pattern of zeros as the upper one. The coefficient matrix has a band structure for the unmodified version of the procedure.

When the basis operators are Hermitian, all elements are purely imaginary. This means that all elements of the main diagonal must be zero. However, if the basis contains non-Hermitian operators, the main diagonal will contain nonzero elements, but they must be real numbers.

3.4 Third step: estimation of the propagator matrix and propagation rules

In principle, the propagator matrix can be obtained by evaluating the matrix exponential according to Eq. (12):

$$\begin{array}{}\text{(24)}& \mathbf{U}=\exp \left(-i\mathbf{L}t\right)=\left(\begin{array}{cccc}{U}_{\mathrm{11}}& U_{\mathrm{12}}& \mathrm{\dots }& U_{\mathrm{1}N}\\ U_{\mathrm{21}}& U_{\mathrm{22}}& \mathrm{\dots }& U_{\mathrm{2}N}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ U_{N\mathrm{1}}& U_{N\mathrm{2}}& \mathrm{\dots }& U_{NN}\end{array}\right).\end{array}$$

L is a pure imaginary matrix for a Hermitian operator basis. In this case, U only contains real elements and is orthogonal. This means that the norms of all the rows and columns of this matrix are unity:

$$\begin{array}{}\text{(25)}& \sum _{i=\mathrm{1}}^N{U}_{ij}^{\mathrm{2}}=\sum _{j=\mathrm{1}}^N{U}_{ij}^{\mathrm{2}}=\mathrm{1}.\end{array}$$

The propagation rules can be obtained from the elements of the transposed propagator matrix (for proof, see Sect. S1.1):

$$\begin{array}{}\text{(26)}& \begin{array}{c}{\widehat{A}}_{\mathrm{1}}\stackrel{\widehat{H}t}{⟶}{\widehat{A}}_{\mathrm{1}}{U}_{\mathrm{11}}+{\widehat{A}}_{\mathrm{2}}{U}_{\mathrm{21}}+\mathrm{\dots }+{\widehat{A}}_N{U}_{N\mathrm{1}},\\ {\widehat{A}}_{\mathrm{2}}\stackrel{\widehat{H}t}{⟶}{\widehat{A}}_{\mathrm{1}}{U}_{\mathrm{12}}+{\widehat{A}}_{\mathrm{2}}{U}_{\mathrm{22}}+\mathrm{\dots }+{\widehat{A}}_N{U}_{N\mathrm{2}},\\ \mathrm{\dots }\end{array}\end{array}$$

However, it is not necessary to recompute the matrix exponential for each new situation. Instead, for low-dimensional subspaces, the generic propagation formulae shown in the next section can be used as a template. Here the elements of the Liouvillian matrix have to be inserted directly, without the need to perform the matrix exponentialization. This was done above for the POF example, where the constant λ resulting from the commutator equations could be used directly as the oscillation frequency.

4 Special cases

The situations in the examples shown below are characterized by different initial states and different Hamiltonians. The Hamiltonians are listed in Appendix B.

A detailed analytical consideration of each example can be found in the Supplement.

4.1 Reduction to a 1-D subspace

In this case, the commutator of the Hamiltonian with the density operator at t=0, i.e., ${\widehat{A}}_{\mathrm{1}}$, is proportional to ${\widehat{A}}_{\mathrm{1}}$ itself,

$$\begin{array}{}\text{(27)}& \left[\widehat{H},{\widehat{A}}_{\mathrm{1}}\right]=\mathit{\lambda }{\widehat{A}}_{\mathrm{1}},\end{array}$$

and only occurs if ${\widehat{A}}_{\mathrm{1}}$ is non-Hermitian, e.g., ${\widehat{A}}_{\mathrm{1}}={\widehat{I}}_{+}\equiv {\widehat{I}}_x+i{\widehat{I}}_y$ as used for the characterization of the complex free induction decay (FID) (Abragam, 1961). Then the coefficient matrix only consists of one scalar λ. The propagator is also a scalar:

$$\begin{array}{}\text{(28)}& {\mathbf{U}}_{\mathrm{1}\mathrm{D}}=e^{-i\mathit{\lambda }t}.\end{array}$$

The corresponding propagation rule is

$$\begin{array}{}\text{(29)}& {\widehat{A}}_{\mathrm{1}}\stackrel{\widehat{H}t}{⟶}{\widehat{A}}_{\mathrm{1}}{e}^{-i\mathit{\lambda }t}.\end{array}$$

Example. A rotating frame, resonance offset Δω, and complex transversal magnetization, which is represented by ${\widehat{I}}_{+}$:

$$\begin{array}{}\text{(30)}& {\widehat{I}}_{+}\stackrel{-\mathrm{\Delta }\mathit{\omega }{\widehat{I}}_{z}t}{⟶}{\widehat{I}}_{+}{e}^{-i\mathrm{\Delta }\mathit{\omega }t}.\end{array}$$

4.2 Case of reduction to a 2-D subspace

The corresponding equations are like Eqs. (3) and (6) and belong to the POF. The Liouvillian matrix and the propagator matrix are the transpositions of the matrices in Eqs. (5) and (6):

$$\begin{array}{}\text{(31)}& \begin{array}{l}{\mathbf{L}}_{\mathrm{2}\mathrm{D}}=\left(\begin{array}{cc}\mathrm{0}& -i\mathit{\lambda }\\ i\mathit{\lambda }& \mathrm{0}\end{array}\right),\\ {\mathbf{U}}_{\mathrm{2}\mathrm{D}}=\exp \left(-i{\mathbf{L}}_{\mathrm{2}\mathrm{D}}t\right)=\left(\begin{array}{cc}\cos \mathit{\lambda }t& -\sin \mathit{\lambda }t\\ \sin \mathit{\lambda }t& \cos \mathit{\lambda }t\end{array}\right).\end{array}\end{array}$$

The propagation formulae are given by Eq. (2), which describes an oscillatory behavior between the initial state and another state described by the commutator of the Hamiltonian, with the operator corresponding to the density operator at the beginning.

In addition to the cases known from numerous POF applications, there are other situations that can be described as time evolution in a 2-D subspace. All of them are well known; they are listed here for the sake of completeness:

• FID of an ensemble of isolated pairs of equal spins (I₁, I₂) after a $\mathit{\pi }/\mathrm{2}$ pulse and homonuclear dipolar interaction within the spin pairs; we observe the transversal magnetization represented by the operator sum ${\widehat{I}}_{\mathrm{1}x}+{\widehat{I}}_{\mathrm{2}x}$:

  $$\begin{array}{}\text{(32)}& \begin{array}{rl}{\widehat{I}}_{\mathrm{1}x}& +{\widehat{I}}_{\mathrm{2}x}\stackrel{{\widehat{H}}_{\mathrm{II}}t}{⟶}\left({\widehat{I}}_{\mathrm{1}x}+{\widehat{I}}_{\mathrm{2}x}\right)\cos {\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{D}_{\mathrm{II}}t\\ & -\mathrm{2}\left({\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}y}+{\widehat{I}}_{\mathrm{1}y}{\widehat{I}}_{\mathrm{2}z}\right)\sin {\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{D}_{\mathrm{II}}t.\end{array}\end{array}$$

• FID of an ensemble of isolated pairs of nonequal spins (I, S) after a $\mathit{\pi }/\mathrm{2}$ pulse in the I channel and heteronuclear dipolar interaction within the spin pairs; we observe the transversal I magnetization represented by the operator ${\widehat{I}}_x$:

  $$\begin{array}{}\text{(33)}& {\widehat{I}}_x\stackrel{{\widehat{H}}_{\mathrm{IS}}t}{⟶}{\widehat{I}}_x\cos D_{\mathrm{IS}}t-\mathrm{2}{\widehat{I}}_y{\widehat{S}}_z\sin D_{\mathrm{IS}}t.\end{array}$$

• FID of an ensemble of spins I=1 (e.g., ²H or ¹⁴N) under quadrupolar interaction; we again follow the transversal magnetization:

  $$\begin{array}{}\text{(34)}& {\widehat{I}}_x\stackrel{{\widehat{H}}_{\mathrm{Q}}t}{⟶}{\widehat{I}}_x\cos {\mathit{\omega }}_{\mathrm{Q}}t+\left({\widehat{I}}_z{\widehat{I}}_y+{\widehat{I}}_y{\widehat{I}}_z\right)\sin {\mathit{\omega }}_{Q}t.\end{array}$$

• Ensemble of pairs of homonuclearly coupled equal spins (I₁, I₂) with spin quantum number $\mathrm{1}/\mathrm{2}$, where initially spin 1 is oriented parallel to B₀ and spin 2 is oriented antiparallel to that. We follow the difference z magnetization, which is represented by ${\widehat{I}}_{\mathrm{1}z}-{\widehat{I}}_{\mathrm{2}z}$:

  $$\begin{array}{}\text{(35)}& \begin{array}{rl}{\widehat{I}}_{\mathrm{1}z}& -{\widehat{I}}_{\mathrm{2}z}\stackrel{{\widehat{H}}_{\mathrm{II}}t}{⟶}\left({\widehat{I}}_{\mathrm{1}z}-{\widehat{I}}_{\mathrm{2}z}\right)\cos D_{\mathrm{II}}t\\ & +\mathrm{2}\left({\widehat{I}}_{\mathrm{1}x}{\widehat{I}}_{\mathrm{2}y}-{\widehat{I}}_{\mathrm{1}y}{\widehat{I}}_{\mathrm{2}x}\right)\sin D_{\mathrm{II}}t.\end{array}\end{array}$$

• Cross-polarization within pairs of antiparallel nonequal spins (I, S): both spins are locked in resonant rf fields with equal nutation frequencies ${\mathit{\omega }}_{\mathrm{1}I}={\mathit{\omega }}_{\mathrm{1}S}\gg D_{\mathrm{IS}}$ (Hartmann–Hahn (HH) condition). The Hamiltonian and the state operators are given in the doubly rotating frame following Hartmann and Hahn (1962), where the z direction is along the rf irradiation. If initially the S spins are oriented parallel to the locking field and the I spins are antiparallel to that, the time evolution can be described by following ${\widehat{S}}_z-{\widehat{I}}_z$:

  $$\begin{array}{}\text{(36)}& \begin{array}{rl}{\widehat{S}}_z& -{\widehat{I}}_z\stackrel{{\widehat{H}}_{\mathrm{HH}}t}{⟶}\left({\widehat{S}}_z-{\widehat{I}}_z\right)\cos D_{\mathrm{IS}}t\\ & +\mathrm{2}\left({\widehat{I}}_x{\widehat{S}}_y-{\widehat{I}}_y{\widehat{S}}_x\right)\sin D_{\mathrm{IS}}t.\end{array}\end{array}$$

  The last equation describes the behavior of the difference between I and S polarizations, not the individual polarizations themselves. The time evolution of the latter requires at least a 3-D approach (see below).

The oscillation takes place in the first three examples between observable transversal magnetization and antiphase states, in the last two examples between longitudinal difference magnetization and zero and double-quantum coherences. These examples show an effect of the dimension reduction: to obtain a 2-D problem, the operators characterizing the states of the spin system have a more complicated structure than in the simple cases above. For example, it would be possible to consider the right-hand side of Eq. (32) to be a linear combination of the four states ${\widehat{I}}_{\mathrm{1}x}$, ${\widehat{I}}_{\mathrm{2}x}$, $\mathrm{2}{\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}y}$, and $\mathrm{2}{\widehat{I}}_{\mathrm{1}y}{\widehat{I}}_{\mathrm{2}z}$ if a more illustrative notation is desired.

4.3 Case of reduction to a 3-D subspace

4.3.1 Generic notation

Here we deal with those cases where the procedure described above reaches the cancellation condition after three commutator equations of the forms

$$\begin{array}{}\text{(37)}& \begin{array}{l}\left[\widehat{H},\widehat{A}\right]=ia\widehat{B},\\ \left[\widehat{H},\widehat{B}\right]=-ia\widehat{A}+ib\widehat{C},\\ \left[\widehat{H},\widehat{C}\right]=-ib\widehat{B},\end{array}\end{array}$$

where a, b∈ℝ. In step 2, we determine the Liouvillian matrix as the transposed coefficient matrix of Eq. (37):

$$\begin{array}{}\text{(38)}& {\mathbf{L}}_{\mathrm{3}\mathrm{D}}=\left(\begin{array}{ccc}\mathrm{0}& -ia& \mathrm{0}\\ ia& \mathrm{0}& -ib\\ \mathrm{0}& ib& \mathrm{0}\end{array}\right).\end{array}$$

From this we determine the matrix of the superpropagator as the matrix exponential corresponding to Eq. (24):

$$\begin{array}{}\text{(39)}& \begin{array}{rl}{\mathbf{U}}_{\mathrm{3}\mathrm{D}}& =\exp \left(-i{\mathbf{L}}_{\mathrm{3}\mathrm{D}}t\right)\\ & ={\displaystyle \frac{\mathrm{1}}{q^{\mathrm{2}}}}\left(\begin{array}{ccc}{b}^{\mathrm{2}}+a^{\mathrm{2}}\cos qt& -aq\sin qt& ab(\mathrm{1}-\cos qt)\\ aq\sin qt& q^{\mathrm{2}}\cos qt& -bq\sin qt\\ ab(\mathrm{1}-\cos qt)& bq\sin qt& a^{\mathrm{2}}+b^{\mathrm{2}}\cos qt\end{array}\right),\end{array}\end{array}$$

with $q^{\mathrm{2}}:=a^{\mathrm{2}}+b^{\mathrm{2}}$. The orthogonality of U_3D, i.e., the validity of Eq. (25), can be verified immediately.

In step (3), according to Eq. (26), the following propagation rules are obtained from the columns of U_3D:

$$\begin{array}{}\text{(40)}& {\displaystyle }\widehat{A}\stackrel{\widehat{H}t}{⟶}\widehat{A}\cdot {\displaystyle \frac{b^{\mathrm{2}}+a^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}+\widehat{B}\cdot {\displaystyle \frac{a}{q}}\sin qt+\widehat{C}\cdot {\displaystyle \frac{ab}{q^{\mathrm{2}}}}(\mathrm{1}-\cos qt),\text{(41)}& {\displaystyle }\widehat{B}\stackrel{\widehat{H}t}{⟶}\widehat{B}\cdot \cos qt+\widehat{C}\cdot {\displaystyle \frac{b}{q}}\sin qt-\widehat{A}\cdot {\displaystyle \frac{a}{q}}\sin qt,\text{(42)}& {\displaystyle }\widehat{C}\stackrel{\widehat{H}t}{⟶}\widehat{C}\cdot {\displaystyle \frac{a^{\mathrm{2}}+b^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}+\widehat{A}\cdot {\displaystyle \frac{ab}{q^{\mathrm{2}}}}(\mathrm{1}-\cos qt)-\widehat{B}\cdot {\displaystyle \frac{b}{q}}\sin qt,\end{array}$$

This can be seen as an extension of the POF to 3-D problems.

If all the basis operators are Hermitian, one of the eigenvalues of L_3D is zero because of detL_3D=0. Therefore, the solution of the Liouville–von Neumann equation may contain a nonzero constant beyond the oscillating terms. The propagator-matrix element ${\left(U_{\mathrm{3}\mathrm{D}}\right)}_{\mathrm{11}}$ contains the time evolution of the initial state $\widehat{A}$. The constant term $b^{\mathrm{2}}/(a^{\mathrm{2}}+b^{\mathrm{2}})$ shows that the oscillations do not take place around zero as in the 2-D case, but around another level. Moreover, its amplitude is reduced to $a^{\mathrm{2}}/(a^{\mathrm{2}}+b^{\mathrm{2}})$, while the frequency increases more the smaller the amplitude is (Fig. 1).

Figure 1 Time evolution of the prefactor of $\widehat{A}$ in Eq. (40) for different ratios $b/a$. The particular case b=0 leads to the 2-D case; therefore, the corresponding curve is a pure oscillation around 0. With increasing b, the average levels of the related oscillations (dashed horizontal lines) increase, whereas the corresponding amplitudes decrease.

There is one important special case: a=b. Here the propagation formulae are simplified to

$$\begin{array}{}\text{(43)}& {\displaystyle }\widehat{A}\stackrel{\widehat{H}t}{⟶}\widehat{A}\cdot {\cos }^{\mathrm{2}}{\displaystyle \frac{qt}{\mathrm{2}}}+\widehat{B}\cdot {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\sin qt+\widehat{C}\cdot {\sin }^{\mathrm{2}}{\displaystyle \frac{qt}{\mathrm{2}}},\text{(44)}& {\displaystyle }\widehat{B}\stackrel{\widehat{H}t}{⟶}\widehat{B}\cdot \cos qt+\widehat{C}\cdot {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\sin qt-\widehat{A}\cdot {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\sin qt,\text{(45)}& {\displaystyle }\widehat{C}\stackrel{\widehat{H}t}{⟶}\widehat{C}\cdot {\cos }^{\mathrm{2}}{\displaystyle \frac{qt}{\mathrm{2}}}+\widehat{A}\cdot {\sin }^{\mathrm{2}}{\displaystyle \frac{qt}{\mathrm{2}}}-\widehat{B}\cdot {\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\sin qt.\end{array}$$

This is applied in the description of cross-polarization and polarization transfer (see the corresponding examples below).

4.3.2 Group 1 of experiments leading to 3-D subspaces: magnetization initially aligned parallel to B₀

• Off-resonance nutation. This involves rf irradiation with strength ω₁ on an ensemble of isolated spins under resonance offset Δω:

  $$\begin{array}{}\text{(46)}& \begin{array}{rl}{\widehat{I}}_z& \stackrel{\left({\widehat{H}}_{\mathrm{Ix}}+{\widehat{H}}_{\mathrm{\Delta }}\right)t}{⟶}{\widehat{I}}_z\cdot {\displaystyle \frac{\mathrm{\Delta }{\mathit{\omega }}^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}\\ & +{\widehat{I}}_y\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}}}{q}}\sin qt+{\widehat{I}}_x\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}}\mathrm{\Delta }\mathit{\omega }}{q^{\mathrm{2}}}}(\mathrm{1}-\cos qt),\end{array}\end{array}$$

  with q² = ${\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{\Delta }{\mathit{\omega }}^{\mathrm{2}}$.

• Nutation and dipolar interaction. The observed spin I, which is heteronuclearly coupled to the spin S, is additionally irradiated with rf of strength ω_1I:

  $$\begin{array}{}\text{(47)}& \begin{array}{rl}{\widehat{I}}_z& \stackrel{\left({\widehat{H}}_{\mathrm{Ix}}+{\widehat{H}}_{\mathrm{IS}}\right)t}{⟶}{\widehat{I}}_z\cdot {\displaystyle \frac{D^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}\\ & +{\widehat{I}}_y\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}}{q}}\sin qt+\mathrm{2}{\widehat{I}}_x{\widehat{S}}_z\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}D}{q^{\mathrm{2}}}}(\mathrm{1}-\cos qt),\end{array}\end{array}$$

  with $q^{\mathrm{2}}={\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+D^{\mathrm{2}}$.

(The case of homonuclear interaction under rf irradiation leads to a 4-D problem; see below.)

The similarity of the two Eqs. (46) and (47) is obvious. For Δω→0 and D_IS→0, respectively, they merge into an equation describing a rotation in the y–z plane. Both equations reflect the well-known fact that a total inversion of the magnetization with a single rectangular pulse is only possible if the offset or the coupling is zero. In addition, the coupling and resonance offset change both the pulse duration ${\mathit{\tau }}_{\mathit{\pi }/\mathrm{2}}$ required to reach maximum y magnetization and the pulse duration τ_π required to reach zero y magnetization to shorter times:

$$\begin{array}{}\text{(48)}& {\mathit{\tau }}_{\mathit{\pi }/\mathrm{2}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}\sqrt{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+C^{\mathrm{2}}}}};{\mathit{\tau }}_{\mathit{\pi }}={\displaystyle \frac{\mathit{\pi }}{\sqrt{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+C^{\mathrm{2}}}}},\end{array}$$

with C=Δω for the off-resonance nutation (Eq. 46) and C=D_IS for the nutation under heteronuclear dipolar interaction.

4.3.3 Group 2 of experiments leading to 3-D subspaces: FID under both rf irradiation and dipolar interaction

• Decoupling experiment: this involves an ensemble of spin pairs {IS}, heteronuclear dipolar interaction between I and S, rf irradiation on the S channel with finite rf power of strength ω_1S, and observations of spin I.

  $$\begin{array}{}\text{(49)}& \begin{array}{rl}{\widehat{I}}_x& \stackrel{\left({\widehat{H}}_{\mathrm{Sx}}+{\widehat{H}}_{\mathrm{IS}}\right)t}{⟶}{\widehat{I}}_x\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}\mathrm{S}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}\\ & -\mathrm{2}{\widehat{I}}_y{\widehat{S}}_z\cdot {\displaystyle \frac{D_{\mathrm{IS}}}{q}}\sin qt\\ & -\mathrm{2}{\widehat{I}}_x{\widehat{S}}_z\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}\mathrm{S}}{D}_{\mathrm{IS}}}{q^{\mathrm{2}}}}(\mathrm{1}-\cos qt),\end{array}\end{array}$$

  with q² = ${\mathit{\omega }}_{\mathrm{1}\mathrm{S}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}$. To obtain the corresponding equation for the J coupling, replace D_IS with −πJ. This equation describes a partial exchange of polarization between x magnetization and two antiphase states.

  The decoupling effect is explained as follows: the observable part of the density operator – the prefactor of ${\widehat{I}}_x$ – contains a constant part ${\mathit{\omega }}_{\mathrm{1}\mathrm{S}}^{\mathrm{2}}/\left({\mathit{\omega }}_{\mathrm{1}\mathrm{S}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\right)$ and an oscillating part with the amplitude $D_{\mathrm{IS}}^{\mathrm{2}}/\left({\mathit{\omega }}_{\mathrm{1}\mathrm{S}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\right)$. Such oscillations were observed for instance in DIPSHIFT experiments (Kurz et al., 2013). Powder-averaging leads to a rather fast decay of the oscillation, which gives a broad line (Pake doublet) after Fourier transformation, while the former gives a δ line. With increasing rf strength ω_1S, the prefactor of the broad peak decreases to zero for infinite rf power, while that of the δ line increases. The constant component is subject to relaxation damping and chemical-shift-induced oscillation on a longer timescale and produces a more or less narrow line.

• On-resonance spin locking and heteronuclear dipolar coupling:

  $$\begin{array}{}\text{(50)}& \begin{array}{rl}{\widehat{I}}_x& \stackrel{\left({\widehat{H}}_{\mathrm{Ix}}+{\widehat{H}}_{\mathrm{IS}}\right)t}{⟶}{\widehat{I}}_x\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}\\ & -\mathrm{2}{\widehat{I}}_y{\widehat{S}}_z\cdot {\displaystyle \frac{D_{\mathrm{IS}}}{q}}\sin qt\\ & +\mathrm{2}{\widehat{I}}_z{\widehat{S}}_z\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}{D}_{\mathrm{IS}}}{q^{\mathrm{2}}}}(\mathrm{1}-\cos qt),\end{array}\end{array}$$

  with q² = ${\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}$. This oscillation frequency is indeed equal to that of the corresponding nutation experiment (Eq. 47).

• On-resonance spin locking and homonuclear dipolar coupling:

  $$\begin{array}{}\text{(51)}& \begin{array}{rl}{\widehat{I}}_{\mathrm{1}\mathrm{x}}& +{\widehat{I}}_{\mathrm{2}\mathrm{x}}\stackrel{\left({\widehat{H}}_{\mathrm{Ix}}+{\widehat{H}}_{\mathrm{II}}\right)t}{⟶}\left({\widehat{I}}_{\mathrm{1}\mathrm{x}}+{\widehat{I}}_{\mathrm{2}\mathrm{x}}\right)\cdot {\displaystyle \frac{\mathrm{4}{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+\frac{\mathrm{9}}{\mathrm{4}}{D}_{\mathrm{II}}^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}\\ & -\mathrm{2}\left({\widehat{I}}_{\mathrm{1}\mathrm{y}}{\widehat{I}}_{\mathrm{2}\mathrm{z}}+{\widehat{I}}_{\mathrm{1}\mathrm{z}}{\widehat{I}}_{\mathrm{2}\mathrm{y}}\right)\cdot {\displaystyle \frac{\mathrm{3}{D}_{\mathrm{II}}}{\mathrm{2}q}}\sin qt\\ & +\mathrm{2}\left({\widehat{I}}_{\mathrm{1}\mathrm{z}}{\widehat{I}}_{\mathrm{2}\mathrm{z}}-{\widehat{I}}_{\mathrm{1}\mathrm{y}}{\widehat{I}}_{\mathrm{2}\mathrm{y}}\right)\cdot {\displaystyle \frac{\mathrm{3}{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}{D}_{\mathrm{II}}}{q^{\mathrm{2}}}}(\mathrm{1}-\cos qt),\end{array}\end{array}$$

  with q² → $\mathrm{4}{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+\frac{\mathrm{9}}{\mathrm{4}}{D}_{\mathrm{IS}}^{\mathrm{2}}$.

Comments on the spin-locking examples (Eqs. 50 and 51):

• Analogous to the other 3-D examples, the oscillation takes place around a level that increases with ω_1I. The latter corresponds to the spin-locked part of the transversal magnetization. At the same time, the amplitude of the oscillation is reduced.

• These oscillations are observed at the beginning of the spin-locking experiments (Krushelnitsky et al., 2018, 2023) and have been described theoretically by Garroway (1979) and McArthur et al. (1969). Due to their orientation dependence, they decay quite rapidly in a powder sample but can be refocused in MAS experiments.

• The propagation Eq. (50) is related to the same Hamiltonian as Eq. (47). As a consequence, the oscillation frequencies are the same. However, the different initial states lead to different subspaces and thus to different propagation formulae.

4.3.4 Group 3 of experiments leading to 3-D subspaces: polarization transfer

The rf field strengths are assumed to be much larger than the corresponding coupling frequencies. This situation is very similar to the polarization transfer treated as a 2-D case. The difference lies in the initial states. Instead of the antiparallel orientation above, here one spin of the pairs is polarized and the other is not. The initial states are now only described by ${\widehat{I}}_{\mathrm{1}\mathrm{z}}$ and ${\widehat{S}}_z$. Since these operators are not elements of the 2-D subspaces of the polarization difference examples above, the motion now takes place in other subspaces, which turn out to be 3-D.

• Equal spins:

  $$\begin{array}{}\text{(52)}& \begin{array}{rl}{\widehat{I}}_{\mathrm{1}\mathrm{z}}& \stackrel{{\widehat{H}}_{\mathrm{II}}t}{⟶}{\widehat{I}}_{\mathrm{1}\mathrm{z}}\cdot {\cos }^{\mathrm{2}}\left({\displaystyle \frac{D_{\mathrm{II}}t}{\mathrm{2}}}\right)\\ & +\left({\widehat{I}}_{\mathrm{1}\mathrm{x}}{\widehat{I}}_{\mathrm{2}\mathrm{y}}+{\widehat{I}}_{\mathrm{1}\mathrm{y}}{\widehat{I}}_{\mathrm{2}\mathrm{x}}\right)\cdot \sin D_{\mathrm{II}}t+{\widehat{I}}_{\mathrm{2}\mathrm{z}}\cdot {\sin }^{\mathrm{2}}\left({\displaystyle \frac{D_{\mathrm{II}}t}{\mathrm{2}}}\right).\end{array}\end{array}$$

• Pair of unequal spins I, S under the Hartmann–Hahn condition in the doubly rotating frame:

  $$\begin{array}{}\text{(53)}& \begin{array}{rl}{\widehat{S}}_z& \stackrel{{\widehat{H}}_{\mathrm{HH}}t}{⟶}{\widehat{S}}_z\cdot {\cos }^{\mathrm{2}}\left({\displaystyle \frac{D_{\mathrm{IS}}t}{\mathrm{2}}}\right)\\ & +\left({\widehat{I}}_x{\widehat{S}}_y+{\widehat{I}}_y{\widehat{S}}_x\right)\cdot \sin D_{\mathrm{IS}}t+{\widehat{I}}_z\cdot {\sin }^{\mathrm{2}}\left({\displaystyle \frac{D_{\mathrm{IS}}t}{\mathrm{2}}}\right).\end{array}\end{array}$$

  Müller et al. (1974) were the first to experimentally demonstrate this oscillatory exchange of polarization during cross-polarization.

• Depolarization of I spins in an ensemble of spin triples $\left\{I,S_{\mathrm{1}},S_{\mathrm{2}}\right\}$ under the Hartmann–Hahn condition in the doubly rotating frame: the coupling frequencies for the I−S₁ and I−S₂ interactions are D₁ and D₂, respectively. The interaction between S₁ and S₂ is assumed to be zero. This can be realized experimentally by irradiating the S spins with a resonance offset which is $\mathrm{1}/\sqrt{\mathrm{2}}$ times the rf strength, known as the Lee–Goldburg condition (Lee and Goldburg, 1965).

  $$\begin{array}{}\text{(54)}& \begin{array}{rl}{\widehat{I}}_z& \stackrel{{\widehat{H}}_{\mathrm{HH}\mathrm{2}}\cdot t}{⟶}{\widehat{I}}_z\cdot {\cos }^{\mathrm{2}}\left({\displaystyle \frac{\sqrt{D_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{2}}^{\mathrm{2}}}}{\mathrm{2}}}t\right)\\ & -{\displaystyle \frac{D_{\mathrm{1}}\left({\widehat{I}}_x{\widehat{S}}_{\mathrm{1}y}-{\widehat{I}}_y{\widehat{S}}_{\mathrm{1}x}\right)+D_{\mathrm{2}}\left({\widehat{I}}_x{\widehat{S}}_{\mathrm{2}y}-{\widehat{I}}_y{\widehat{S}}_{\mathrm{2}x}\right)}{\sqrt{D_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{2}}^{\mathrm{2}}}}}\sin \sqrt{D_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{2}}^{\mathrm{2}}}t\\ & +{\displaystyle \frac{D_{\mathrm{1}}^{\mathrm{2}}{\widehat{S}}_{\mathrm{1}z}+D_{\mathrm{2}}^{\mathrm{2}}{\widehat{S}}_{\mathrm{2}z}-\mathrm{4}{D}_{\mathrm{1}}{D}_{\mathrm{2}}{\widehat{I}}_z\left({\widehat{S}}_{\mathrm{1}x}{\widehat{S}}_{\mathrm{2}x}+{\widehat{S}}_{\mathrm{1}y}{\widehat{S}}_{\mathrm{2}y}\right)}{D_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{2}}^{\mathrm{2}}}}{\sin }^{\mathrm{2}}\left({\displaystyle \frac{\sqrt{D_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{2}}^{\mathrm{2}}}}{\mathrm{2}}}t\right)\end{array}\end{array}$$

  Comments on the third item: the oscillation frequency is the geometric sum of the individual frequencies. The oscillation takes place between the initial state and a mixture of observable and unobservable states.

Note the first two propagation rules of the third group: in addition to the orthogonality relations (Eq. 25), the linear sum of the prefactors of the first and third terms is 1. Both terms represent z magnetization. This condition reflects the fact that the sum of the z polarizations is constant for cross-polarization and polarization transfer. This is supported by the fact that ${\widehat{I}}_{\mathrm{1}\mathrm{z}}+{\widehat{I}}_{\mathrm{2}\mathrm{z}}$ commutes with ${\widehat{H}}_{\mathrm{II}}$ and ${\widehat{I}}_z+{\widehat{S}}_z$ commutes with ${\widehat{H}}_{\mathrm{HH}}$. Figure 2 shows the time evolution of the three prefactors in these propagation rules.

Figure 2 Time evolution of the three prefactors in the propagation rules shown in Eqs. (52) and (53).

4.3.5 Group 4 of 3-D examples: cross-polarization, finite rf power, and possible deviation from the HH condition

• Considering the polarization difference in the rotating frame,

  $$\begin{array}{}\text{(55)}& \begin{array}{rl}{\widehat{S}}_x& -{\widehat{I}}_x\stackrel{\left({\widehat{H}}_{\mathrm{Ix}}+{\widehat{H}}_{\mathrm{Sx}}+{\widehat{H}}_{\mathrm{IS}}\right)t}{⟶}\left({\widehat{S}}_x-{\widehat{I}}_x\right)\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{\Delta }}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos q_{\mathrm{\Delta }}t}{q_{\mathrm{\Delta }}^{\mathrm{2}}}}\\ & -\mathrm{2}\left({\widehat{I}}_z{\widehat{S}}_y-{\widehat{I}}_y{\widehat{S}}_z\right)\cdot {\displaystyle \frac{D_{\mathrm{IS}}}{q_{\mathrm{\Delta }}}}\sin q_{\mathrm{\Delta }}t\\ & +\mathrm{2}\left({\widehat{I}}_z{\widehat{S}}_z+{\widehat{I}}_y{\widehat{S}}_y\right)\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{\Delta }}{D}_{\mathrm{IS}}}{q_{\mathrm{\Delta }}^{\mathrm{2}}}}{\sin }^{\mathrm{2}}{q}_{\mathrm{\Delta }}t,\end{array}\end{array}$$

  with ω_Δ = ω_1S−ω_1I and $q_{\mathrm{\Delta }}^{\mathrm{2}}$ = $D_{\mathrm{IS}}^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{\Delta }}^{\mathrm{2}}$, there are the following comments:

  • -

    The CP oscillation frequency is no longer D_IS as calculated for infinite rf power (see above), but $\sqrt{D_{\mathrm{IS}}^{\mathrm{2}}+{\left({\mathit{\omega }}_{\mathrm{1}\mathrm{I}}-{\mathit{\omega }}_{\mathrm{1}\mathrm{S}}\right)}^{\mathrm{2}}}$, which increases with the difference in the two rf field strengths.

  • -

    Moreover, only the relative part $D_{\mathrm{IS}}^{\mathrm{2}}/q_{\mathrm{\Delta }}^{\mathrm{2}}$ of the total magnetization participates in the oscillation. That is, the greater the deviation from the Hartmann–Hahn condition, the lower the maximum transmitted polarization.

• Consider the polarization sum:

  $$\begin{array}{}\text{(56)}& \begin{array}{rl}{\widehat{S}}_x& +{\widehat{I}}_x\stackrel{\left({\widehat{H}}_{\mathrm{Ix}}+{\widehat{H}}_{\mathrm{Sx}}+{\widehat{H}}_{\mathrm{IS}}\right)t}{⟶}\left({\widehat{S}}_x+{\widehat{I}}_x\right)\cdot {\displaystyle \frac{\mathrm{4}{\mathit{\omega }}_{\mathrm{ø}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos q_{\mathrm{ø}}t}{q_{\mathrm{ø}}^{\mathrm{2}}}}\\ & -\mathrm{2}\left({\widehat{I}}_z{\widehat{S}}_y+{\widehat{I}}_y{\widehat{S}}_z\right)\cdot {\displaystyle \frac{D_{\mathrm{IS}}}{q_{\mathrm{ø}}}}\sin q_{\mathrm{ø}}t\\ & +\mathrm{2}\left({\widehat{I}}_z{\widehat{S}}_z-{\widehat{I}}_y{\widehat{S}}_y\right)\cdot {\displaystyle \frac{\mathrm{2}{\mathit{\omega }}_{\mathrm{ø}}{D}_{\mathrm{IS}}}{q_{\mathrm{ø}}^{\mathrm{2}}}}{\sin }^{\mathrm{2}}{q}_{\mathrm{ø}}t,\end{array}\end{array}$$

  with ω_ø = $\left({\mathit{\omega }}_{\mathrm{1}\mathrm{S}}+{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}\right)/\mathrm{2}$ and $q_{\mathrm{ø}}^{\mathrm{2}}$ = $D_{\mathrm{IS}}^{\mathrm{2}}+\mathrm{4}{\mathit{\omega }}_{\mathrm{ø}}^{\mathrm{2}}$.

In contrast to the relations obtained for infinite ω₁, the sum of both polarizations oscillates. The amplitude decreases with increasing rf power. This phenomenon is analogous to what happens with spin locking and decoupling (see the corresponding examples above).

4.4 Case of reduction to a 4-D subspace

4.4.1 Generic notation

This group of situations can be described by commutator relations of the forms

$$\begin{array}{}\text{(57)}& \begin{array}{rl}\left[\widehat{H},\widehat{A}\right]& =ia\widehat{B},\\ \left[\widehat{H},\widehat{B}\right]& =-ia\widehat{A}+ib\widehat{C},\\ \left[\widehat{H},\widehat{C}\right]& =-ib\widehat{B}\mp ia\widehat{D},\\ \left[\widehat{H},\widehat{D}\right]& =\pm ia\widehat{C},\end{array}\end{array}$$

where $\widehat{A}$, $\widehat{B}$, $\widehat{C}$, and $\widehat{D}$ are pairwise orthogonal operators with equal norms. According to step 2, we obtain the Liouvillian matrix from these rules as a transposed coefficient matrix

$$\begin{array}{}\text{(58)}& {\mathbf{L}}_{\mathrm{4}\mathrm{D}}=\left(\begin{array}{cccc}\mathrm{0}& -ia& \mathrm{0}& \mathrm{0}\\ ia& \mathrm{0}& -ib& \mathrm{0}\\ \mathrm{0}& ib& \mathrm{0}& \pm ia\\ \mathrm{0}& \mathrm{0}& \mp ia& \mathrm{0}\end{array}\right)\end{array}$$

and the superpropagator matrix as

$$\begin{array}{}\text{(59)}& \begin{array}{rl}& {\mathbf{U}}_{\mathrm{4}\mathrm{D}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}W}}\\ & \left(\begin{array}{cccc}{q}_{\mathrm{2}}\cos q_{\mathrm{1}}t+q_{\mathrm{1}}\cos q_{\mathrm{2}}t& -a(\sin q_{\mathrm{1}}t+\sin q_{\mathrm{2}}t)& a(\cos q_{\mathrm{1}}t-\cos q_{\mathrm{2}}t)& \mp q_{\mathrm{2}}\sin q_{\mathrm{1}}t\pm q_{\mathrm{1}}\sin q_{\mathrm{2}}t\\ a(\sin q_{\mathrm{1}}t+\sin q_{\mathrm{2}}t)& q_{\mathrm{1}}\cos q_{\mathrm{1}}t+q_{\mathrm{2}}\cos q_{\mathrm{2}}t& q_{\mathrm{2}}\sin q_{\mathrm{2}}t-q_{\mathrm{1}}\sin q_{\mathrm{1}}t& \pm a(\cos q_{\mathrm{2}}t-\cos q_{\mathrm{1}}t)\\ a(\cos q_{\mathrm{2}}t-\cos q_{\mathrm{1}}t)& q_{\mathrm{1}}\sin q_{\mathrm{1}}t-q_{\mathrm{2}}\sin q_{\mathrm{2}}t& q_{\mathrm{1}}\cos q_{\mathrm{1}}t+q_{\mathrm{2}}\cos q_{\mathrm{2}}t& \mp a(\sin q_{\mathrm{1}}t+\sin q_{\mathrm{2}}t)\\ \pm q_{\mathrm{2}}\sin q_{\mathrm{1}}t\mp q_{\mathrm{1}}\sin q_{\mathrm{2}}t& \pm a(\cos q_{\mathrm{2}}t-\cos q_{\mathrm{1}}t)& \pm a(\sin q_{\mathrm{1}}t+\sin q_{\mathrm{2}}t)& q_{\mathrm{2}}\cos q_{\mathrm{1}}t+q_{\mathrm{1}}\cos q_{\mathrm{2}}t\end{array}\right),\end{array}\end{array}$$

with $W^{\mathrm{2}}:=(a^{\mathrm{2}}+b^{\mathrm{2}})/\mathrm{4}$ and $q_{\mathrm{1};\mathrm{2}}:=W\mp b/\mathrm{2}$. The propagation rule for the case where $\widehat{A}$ was the initial state can be read from the first column of U_4D in Eq. (59):

$$\begin{array}{}\text{(60)}& \begin{array}{rl}\widehat{A}\stackrel{\widehat{H}\cdot t}{⟶}& \widehat{A}\cdot {\displaystyle \frac{q_{\mathrm{2}}\cos q_{\mathrm{1}}t+q_{\mathrm{1}}\cos q_{\mathrm{2}}t}{\mathrm{2}W}}\\ & +\widehat{B}\cdot {\displaystyle \frac{a(\sin q_{\mathrm{1}}t+\sin q_{\mathrm{2}}t)}{\mathrm{2}W}}\\ & +\widehat{C}\cdot {\displaystyle \frac{a(\cos q_{\mathrm{2}}t-\cos q_{\mathrm{1}}t)}{\mathrm{2}W}}\\ & \pm \widehat{D}\cdot {\displaystyle \frac{q_{\mathrm{2}}\sin q_{\mathrm{1}}t-q_{\mathrm{1}}\sin q_{\mathrm{2}}t}{\mathrm{2}W}}.\end{array}\end{array}$$

In some cases, another form of this propagation formula may be appropriate for use; this is obtained from Eq. (60) by applying some trigonometric rules:

$$\begin{array}{}\text{(61)}& {\displaystyle }\begin{array}{rl}\widehat{A}\stackrel{\widehat{H}\cdot t}{⟶}& \widehat{A}\cdot \left(\cos \left(Wt\right)\cos {\displaystyle \frac{bt}{\mathrm{2}}}+{\displaystyle \frac{b}{\mathrm{2}W}}\sin \left(Wt\right)\sin {\displaystyle \frac{bt}{\mathrm{2}}}\right)\\ & +\widehat{B}\cdot {\displaystyle \frac{a}{W}}\sin \left(Wt\right)\cos {\displaystyle \frac{bt}{\mathrm{2}}}\end{array}\text{(62)}& {\displaystyle }\begin{array}{rl}& -\widehat{C}\cdot {\displaystyle \frac{a}{W}}\sin \left(Wt\right)\sin {\displaystyle \frac{bt}{\mathrm{2}}}\\ & \mp \widehat{D}\cdot \left(\cos \left(Wt\right)\sin {\displaystyle \frac{bt}{\mathrm{2}}}+{\displaystyle \frac{b}{\mathrm{2}W}}\sin \left(Wt\right)\cos {\displaystyle \frac{bt}{\mathrm{2}}}\right).\end{array}\end{array}$$

4.4.2 Examples

• AB spin system with J coupling and distance of the two lines $\mathrm{\Delta }\mathit{\nu }:=\mathrm{\Delta }\mathit{\omega }/\left(\mathrm{2}\mathit{\pi }\right)$; the spectrometer frequency is assumed to be set at the midpoint between the two resonances:

  $$\begin{array}{}\text{(63)}& \begin{array}{rl}& {\widehat{I}}_{\mathrm{1}x}+{\widehat{I}}_{\mathrm{2}x}\stackrel{\mathrm{2}\mathit{\pi }J{\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}z}t}{⟶}\\ & \left({\widehat{I}}_{\mathrm{1}x}+{\widehat{I}}_{\mathrm{2}x}\right)\cdot {\displaystyle \frac{q_{\mathrm{2}}\cos q_{\mathrm{1}}t+q_{\mathrm{1}}\cos q_{\mathrm{2}}t}{W}}\\ & +\left({\widehat{I}}_{\mathrm{1}y}-{\widehat{I}}_{\mathrm{2}y}\right)\cdot {\displaystyle \frac{\mathrm{\Delta }\mathit{\omega }}{\mathrm{2}}}{\displaystyle \frac{\sin q_{\mathrm{1}}t+\sin q_{\mathrm{2}}t}{W}}\\ & +\mathrm{2}\left({\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}x}-{\widehat{I}}_{\mathrm{1}x}{\widehat{I}}_{\mathrm{2}z}\right)\cdot {\displaystyle \frac{\mathrm{\Delta }\mathit{\omega }}{\mathrm{2}}}{\displaystyle \frac{\cos q_{\mathrm{1}}t-\cos q_{\mathrm{2}}t}{W}}\\ & -\mathrm{2}\left({\widehat{I}}_{\mathrm{1}y}{\widehat{I}}_{\mathrm{2}z}+{\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}y}\right)\cdot {\displaystyle \frac{q_{\mathrm{1}}\sin q_{\mathrm{2}}t-q_{\mathrm{2}}\sin q_{\mathrm{1}}t}{W}},\end{array}\end{array}$$

  where $W=\mathrm{2}\mathit{\pi }\sqrt{J^{\mathrm{2}}+\mathrm{\Delta }{\mathit{\nu }}^{\mathrm{2}}}$ and $q_{\mathrm{1};\mathrm{2}}=\mathit{\pi }\left(\sqrt{J^{\mathrm{2}}+\mathrm{\Delta }{\mathit{\nu }}^{\mathrm{2}}}\mp J\right)$. The cosine terms contain two frequencies that give four line positions ±q₁ and ±q₂ symmetrically around zero after complex Fourier transformation. The intensities are given by the prefactors of the corresponding trigonometric functions. The lower-frequency oscillation has the larger prefactor, which reflects the roof effect. See, e.g., Abragam (1961), Chap. XI, Sect. B. In this book, positions and intensities are calculated from transition frequencies and probabilities for the transitions between the levels. Figure 3 shows the relationship between q₁, q₂, and the intensities and positions of the four lines of an AB spin system.

  For the case Δν=0, the Hamiltonian commutes with the initial state, with the consequence that the density operator remains constant. The Fourier transform of this is just a resonance at zero frequency.

• rf irradiation onto homonuclearly coupled spins $\mathit{1}/\mathit{2}$ which are initially in equilibrium (nutation):

  $$\begin{array}{}\text{(64)}& \begin{array}{rl}& {\widehat{I}}_{\mathrm{1}\mathrm{z}}+{\widehat{I}}_{\mathrm{2}\mathrm{z}}\stackrel{\left({\widehat{H}}_{\mathrm{II}}+{\widehat{H}}_{\mathrm{Ix}}\right)t}{⟶}\left({\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z}\right)\\ & \cdot \left(\cos Wt\cos {\displaystyle \frac{\mathrm{3}{D}_{\mathrm{II}}}{\mathrm{4}}}t+{\displaystyle \frac{\mathrm{3}{D}_{II}}{\mathrm{4}W}}\sin Wt\sin {\displaystyle \frac{\mathrm{3}{D}_{\mathrm{II}}}{\mathrm{4}}}t\right)\\ & +\left({\widehat{I}}_{\mathrm{1}y}+{\widehat{I}}_{\mathrm{2}y}\right)\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}}{W}}\sin Wt\cos {\displaystyle \frac{\mathrm{3}{D}_{\mathrm{II}}}{\mathrm{4}}}t\\ & -\mathrm{2}\left({\widehat{I}}_{\mathrm{1}x}{\widehat{I}}_{\mathrm{2}z}+{\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}x}\right)\cdot {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}}{W}}\sin Wt\cos {\displaystyle \frac{\mathrm{3}{D}_{\mathrm{II}}}{\mathrm{4}}}t\\ & -\mathrm{2}\left({\widehat{I}}_{\mathrm{1}y}{\widehat{I}}_{\mathrm{2}x}+{\widehat{I}}_{\mathrm{1}x}{\widehat{I}}_{\mathrm{2}y}\right)\cdot \left(\cos Wt\sin {\displaystyle \frac{\mathrm{3}{D}_{\mathrm{II}}}{\mathrm{4}}}t-{\displaystyle \frac{\mathrm{3}{D}_{\mathrm{II}}}{\mathrm{4}W}}\sin Wt\cos {\displaystyle \frac{\mathrm{3}{D}_{\mathrm{II}}}{\mathrm{4}}}t\right),\end{array}\end{array}$$

  where W = $\sqrt{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+\frac{\mathrm{9}}{\mathrm{16}}{D}_{\mathrm{II}}^{\mathrm{2}}}$.

  This propagation formula describes the effect of a limited-power rf pulse on the equilibrium magnetization. The time evolution of the prefactor of ${\widehat{I}}_{\mathrm{1}y}+{\widehat{I}}_{\mathrm{2}y}$ is shown in Fig. 4.

  Comments on Eq. (64):

  • -

    As the coupling frequency increases, so does the nutation frequency W. However, this oscillation is modulated by half of the dipolar frequency (Fig. 4).

  • -

    As a consequence, the π and $\mathit{\pi }/\mathrm{2}$ conditions for achieving maximum and zero y magnetization are modified with respect to the coupling-free case. Similar to Eq. (48), this results in

    $$\begin{array}{}\text{(65)}& \begin{array}{rl}& {\mathit{\tau }}_{\mathit{\pi }/\mathrm{2}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}\sqrt{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+\frac{\mathrm{9}}{\mathrm{16}}{D}_{\mathrm{II}}^{\mathrm{2}}}}},\\ & {\mathit{\tau }}_{\mathit{\pi }}={\displaystyle \frac{\mathit{\pi }}{\sqrt{{\mathit{\omega }}_{\mathrm{1}\mathrm{I}}^{\mathrm{2}}+\frac{\mathrm{9}}{\mathrm{16}}{D}_{\mathrm{II}}^{\mathrm{2}}}}}.\end{array}\end{array}$$

• Nutation under quadrupolar interaction, spin 1:

  $$\begin{array}{}\text{(66)}& \begin{array}{rl}& {\widehat{I}}_z\stackrel{\left({\widehat{H}}_Q+{\widehat{H}}_{\mathrm{Ix}}\right)t}{⟶}\\ & {\widehat{I}}_z\cdot \left(\cos {\displaystyle \frac{Wt}{\mathrm{2}}}\cos {\displaystyle \frac{{\mathit{\omega }}_{\mathrm{Q}}t}{\mathrm{2}}}+{\displaystyle \frac{{\mathit{\omega }}_{Q}t}{W}}\sin {\displaystyle \frac{Wt}{\mathrm{2}}}\sin {\displaystyle \frac{{\mathit{\omega }}_{Q}t}{\mathrm{2}}}\right)\\ & +{\widehat{I}}_y\cdot {\displaystyle \frac{\mathrm{2}{\mathit{\omega }}_{\mathrm{1}}}{W}}\sin {\displaystyle \frac{Wt}{\mathrm{2}}}\sin {\displaystyle \frac{{\mathit{\omega }}_{Q}t}{\mathrm{2}}}-\left({\widehat{I}}_x{\widehat{I}}_z+{\widehat{I}}_z{\widehat{I}}_x\right)\\ & \cdot {\displaystyle \frac{\mathrm{2}{\mathit{\omega }}_{\mathrm{1}}}{W}}\sin {\displaystyle \frac{Wt}{\mathrm{2}}}\sin {\displaystyle \frac{{\mathit{\omega }}_{Q}t}{\mathrm{2}}}+\left({\widehat{I}}_y{\widehat{I}}_x+{\widehat{I}}_x{\widehat{I}}_y\right)\\ & \cdot \left(\cos {\displaystyle \frac{Wt}{\mathrm{2}}}\sin {\displaystyle \frac{{\mathit{\omega }}_{Q}t}{\mathrm{2}}}+{\displaystyle \frac{{\mathit{\omega }}_{Q}t}{W}}\sin {\displaystyle \frac{Wt}{\mathrm{2}}}\cos {\displaystyle \frac{{\mathit{\omega }}_{Q}t}{\mathrm{2}}}\right),\end{array}\end{array}$$

  with W = $\sqrt{\mathrm{4}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\omega }}_Q^{\mathrm{2}}}$. This is consistent with the findings of Bloom et al. (1980), Barbara et al. (1986), and Vega and Luz (1987). Again, the result is that nutation occurs more quickly than ω_1I if there is an additional interaction.

Figure 3 Connection between the two oscillation frequencies, the positions, and the intensities of lines in the spectrum of an AB spin system (example: Δν=3J).

Figure 4 Nutation curves for different ratios of rf field strength and coupling strength. For homonuclear dipolar coupling, $D_{\mathrm{II}}/{\mathit{\omega }}_{\mathrm{1}}$ = (a) 0, (b) 0.15, (c) 0.3, (d) 0.45, (e) 0.6, and (f) 0.75. For quadrupolar coupling, ${\mathit{\omega }}_Q/{\mathit{\omega }}_{\mathrm{1}}$ = (a) 0, (b) 0.2, (c) 0.4, (d) 0.6, (e) 0.8, and (f) 1.

4.5 Example of a 5-D subspace

We consider cross-polarization under a finite rf power, assuming that the Hartmann–Hahn condition is satisfied: ${\mathit{\omega }}_{\mathrm{1}\mathrm{I}}={\mathit{\omega }}_{\mathrm{1}\mathrm{S}}=:{\mathit{\omega }}_{\mathrm{1}}$. Other than in the corresponding 3-D examples shown above, the initial state consists of a transversally polarized S spin and a depolarized I spin, i.e., $\widehat{\mathit{\rho }}\left(\mathrm{0}\right)={\widehat{S}}_x$:

$$\begin{array}{}\text{(67)}& \begin{array}{rl}{\widehat{S}}_x& \stackrel{({\widehat{H}}_{\mathrm{Ix}}+{\widehat{H}}_{\mathrm{Sx}}+{\widehat{H}}_{\mathrm{IS}})\cdot t}{⟶}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}[\left(\cos D_{\mathrm{IS}}t+{\displaystyle \frac{\mathrm{4}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}\right)\cdot {\widehat{S}}_x\\ & -\left({\displaystyle \frac{D_{\mathrm{IS}}}{q}}\sin qt+\sin D_{\mathrm{IS}}t\right)\cdot \mathrm{2}{\widehat{I}}_z{\widehat{S}}_y\\ & +\mathrm{4}{\displaystyle \frac{D_{\mathrm{IS}}{\mathit{\omega }}_{\mathrm{1}}}{q^{\mathrm{2}}}}(\mathrm{1}-\cos qt)\left({\widehat{I}}_z{\widehat{S}}_z-{\widehat{I}}_y{\widehat{S}}_y\right)\\ & +\left({\displaystyle \frac{D_{\mathrm{IS}}}{q}}\sin qt-\sin D_{\mathrm{IS}}t\right)\cdot \mathrm{2}{\widehat{I}}_y{\widehat{S}}_z\\ & +\left(-\cos D_{\mathrm{IS}}t+{\displaystyle \frac{\mathrm{4}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}\right)\cdot {\widehat{I}}_x].\end{array}\end{array}$$

Equation (67) describes the evolution of the magnetization of the initially polarized spin. That of the other spin evolves as given in Eq. (68):

$$\begin{array}{}\text{(68)}& \begin{array}{rl}{\widehat{I}}_x& \stackrel{({\widehat{H}}_{\mathrm{Ix}}+{\widehat{H}}_{\mathrm{Sx}}+{\widehat{H}}_{\mathrm{IS}})\cdot t}{⟶}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}[{\widehat{S}}_x\cdot \left(-\cos D_{\mathrm{IS}}t+{\displaystyle \frac{\mathrm{4}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}\right)\\ & -\left({\displaystyle \frac{D_{\mathrm{IS}}}{q}}\sin qt-\sin D_{\mathrm{IS}}t\right)\cdot \mathrm{2}{\widehat{I}}_z{\widehat{S}}_y\\ & +\mathrm{4}{\displaystyle \frac{D_{\mathrm{IS}}{\mathit{\omega }}_{\mathrm{1}}}{q^{\mathrm{2}}}}(\mathrm{1}-\cos qt)\left({\widehat{I}}_z{\widehat{S}}_z-{\widehat{I}}_y{\widehat{S}}_y\right)\\ & -\left({\displaystyle \frac{D_{\mathrm{IS}}}{q}}\sin qt+\sin D_{\mathrm{IS}}t\right)\cdot \mathrm{2}{\widehat{I}}_y{\widehat{S}}_z\\ & +\left(\cos D_{\mathrm{IS}}t+{\displaystyle \frac{\mathrm{4}{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos qt}{q^{\mathrm{2}}}}\right)\cdot {\widehat{I}}_x].\end{array}\end{array}$$

4.6 Examples of a 6-D subspace

• In the Hartmann–Hahn cross-polarization experiment with finite rf power and deviation from the Hartmann–Hahn condition, the initial state consists of a transversally polarized S spin and a depolarized I spin, i.e., $\widehat{\mathit{\rho }}\left(\mathrm{0}\right)={\widehat{S}}_x$. Using the variables defined in Sect. 4.3.5, we obtain

  $$\begin{array}{}\text{(69)}& \begin{array}{rl}& {\widehat{S}}_x\stackrel{({\widehat{H}}_{\mathrm{Ix}}+{\widehat{H}}_{\mathrm{Sx}}+{\widehat{H}}_{\mathrm{IS}})\cdot t}{⟶}\\ & {\widehat{S}}_x\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{{\mathit{\omega }}_{\mathrm{\Delta }}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos q_{\mathrm{\Delta }}t}{q_{\mathrm{\Delta }}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{4}{\mathit{\omega }}_{\mathrm{ø}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos q_{\mathrm{ø}}t}{q_{\mathrm{ø}}^{\mathrm{2}}}}\right)\\ & -{\widehat{I}}_z{\widehat{S}}_y\cdot D_{\mathrm{IS}}\left({\displaystyle \frac{\sin q_{\mathrm{\Delta }}t}{q_{\mathrm{\Delta }}}}+{\displaystyle \frac{\sin q_{\mathrm{ø}}t}{q_{\mathrm{ø}}}}\right)\\ & +\mathrm{2}{\widehat{I}}_z{\widehat{S}}_z{D}_{\mathrm{IS}}\left[{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{\Delta }}}{\mathrm{2}{q}_{\mathrm{\Delta }}}}(\mathrm{1}-\cos q_{\mathrm{\Delta }}t)+{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{ø}}}{q_{\mathrm{ø}}}}(\mathrm{1}-\cos q_{\mathrm{ø}}t)\right]\\ & +\mathrm{2}{\widehat{I}}_y{\widehat{S}}_y{D}_{\mathrm{IS}}\left[{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{\Delta }}}{\mathrm{2}{q}_{\mathrm{\Delta }}}}(\mathrm{1}-\cos q_{\mathrm{\Delta }}t)-{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{ø}}}{q_{\mathrm{ø}}}}(\mathrm{1}-\cos q_{\mathrm{ø}}t)\right]\\ & +{\widehat{I}}_y{\widehat{S}}_z\cdot D_{\mathrm{IS}}\left({\displaystyle \frac{\sin q_{\mathrm{\Delta }}t}{q_{\mathrm{\Delta }}}}-{\displaystyle \frac{\sin q_{\mathrm{ø}}t}{q_{\mathrm{ø}}}}\right)\\ & +{\widehat{I}}_x\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{4}{\mathit{\omega }}_{\mathrm{ø}}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos q_{\mathrm{ø}}t}{q_{\mathrm{ø}}^{\mathrm{2}}}}-{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{\Delta }}^{\mathrm{2}}+D_{\mathrm{IS}}^{\mathrm{2}}\cos q_{\mathrm{\Delta }}t}{q_{\mathrm{\Delta }}^{\mathrm{2}}}}\right).\end{array}\end{array}$$

  Figure 5 shows the buildup curve for the I magnetization for the case where it was initially unpolarized and S was polarized. For zero deviation from the Hartmann–Hahn condition, the curve still looks similar to the squared sine derived for infinite rf power (Eq. 53). Small deformations result from the fact that the rf power is finite. An increasing deviation from the Hartmann–Hahn condition leads to a strong loss of the polarization-transfer efficiency, in addition to further deviation from the ideal curve.

• Hartmann–Hahn cross-polarization of an I spin from two S spins. The Hamiltonian is the same as in Sect. 4.3.4, and we consider the problem in the doubly rotating frame. In this example, however, we want to follow the time evolution of all three spins individually which could not be separated in the 3-D example. The propagation rules for all three spins can be found in the Supplement. Here is the propagation rule for the case where the system starts with polarized S spins, i.e., $\widehat{A}={\widehat{S}}_{\mathrm{1}\mathrm{z}}+{\widehat{S}}_{\mathrm{2}\mathrm{z}}$ and $\widehat{H}={\widehat{H}}_{\mathrm{HH}\mathrm{2}}$:

  $$\begin{array}{}\text{(70)}& \begin{array}{rl}& {\widehat{S}}_{\mathrm{1}\mathrm{z}}+{\widehat{S}}_{\mathrm{2}\mathrm{z}}\stackrel{{\widehat{H}}_{\mathrm{HH}\mathrm{2}}\cdot t}{⟶}{\widehat{I}}_z\cdot {\sin }^{\mathrm{2}}{\displaystyle \frac{qt}{\mathrm{2}}}\\ & +\left({\widehat{S}}_{\mathrm{1}\mathrm{z}}+{\widehat{S}}_{\mathrm{2}\mathrm{z}}\right)\cdot {\displaystyle \frac{\mathrm{1}+{\cos }^{\mathrm{2}}\frac{qt}{\mathrm{2}}}{\mathrm{2}}}\\ & +\left({\widehat{S}}_{\mathrm{1}\mathrm{z}}-{\widehat{S}}_{\mathrm{2}\mathrm{z}}\right)\cdot {\displaystyle \frac{D_{\mathrm{2}}^{\mathrm{2}}-D_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}{q}^{\mathrm{2}}}}{\sin }^{\mathrm{2}}{\displaystyle \frac{qt}{\mathrm{2}}}\\ & +\left({\widehat{I}}_x{\widehat{S}}_{\mathrm{1}\mathrm{y}}-{\widehat{I}}_y{\widehat{S}}_{\mathrm{1}\mathrm{x}}\right){\displaystyle \frac{D_{\mathrm{1}}}{q}}\sin qt\\ & +\left({\widehat{I}}_x{\widehat{S}}_{\mathrm{2}\mathrm{y}}-{\widehat{I}}_y{\widehat{S}}_{\mathrm{2}\mathrm{x}}\right){\displaystyle \frac{D_{\mathrm{2}}}{q}}\sin qt\\ & -{\widehat{I}}_z\left({\widehat{S}}_{\mathrm{1}\mathrm{x}}{\widehat{S}}_{\mathrm{2}\mathrm{x}}-{\widehat{S}}_{\mathrm{1}\mathrm{y}}{\widehat{S}}_{\mathrm{2}\mathrm{y}}\right)\cdot \mathrm{4}{\displaystyle \frac{D_{\mathrm{1}}{D}_{\mathrm{2}}}{q^{\mathrm{2}}}}{\sin }^{\mathrm{2}}{\displaystyle \frac{qt}{\mathrm{2}}},\end{array}\end{array}$$

  with q = $\sqrt{D_{\mathrm{1}}^{\mathrm{2}}+D_{\mathrm{2}}^{\mathrm{2}}}$. This propagation rule describes the oscillatory polarization exchange between the spin polarizations and three unobservable states on the one hand and the oscillatory polarization transfer between the spins on the other hand, as may happen, for example, in ¹³CH₂ or ¹⁵NH₂ groups.

Figure 5 Examples of cross-polarization buildup curves for the I spins in an ensemble of spin pairs I, S corresponding to Eq. (69) with ω_ø=3D_IS. The deviation ω_Δ from the Hartmann–Hahn condition is varied.

5 Outlook: time-dependent Hamiltonians and Liouvillians

In order to obtain the propagator matrices and the propagation formulae of the respective Liouvillians, it was assumed in the above sections that the interactions and thus also the Hamiltonian and Liouvillian are time-invariant. However, the method can be extended to situations where the interaction constants vary with time, fluctuating due to thermal motion or periodically due to sample spinning. The commutator equations are also valid; i.e., steps (1) and (2) of the given procedure can be performed to obtain the relevant subspace and the Liouvillian matrix associated with that subspace. In other words, the Liouvillian matrices obtained for the above examples (see the Supplement) can also be used in the time-dependent situations. The Liouville–von Neumann equation now belongs to a system of linear differential equations but with time-dependent coefficients. There is no general scheme for their integration. The propagator matrix is no longer the matrix exponential of −iLt.

However, in some cases solutions are possible in the following ways: (i) use of the time-averaged Liouvillian for an exact solution when L depends linearly on a single time-dependent parameter and (ii) use of the Shirley method, based on the Floquet theorem, also for an analytical solution.

The first way can be justified as follows: for the Liouville–von Neumann equation in matrix form, one can try to find an effective Liouvillian matrix L_eff(t), which is defined as that matrix which is constant in the interval (0,t) and which has the same effect as the actual L, i.e., $\mathbf{U}\left(t\right)=\exp (-i{\mathbf{L}}_{\mathrm{eff}}t)$ in matrix notation. This can be done using the Magnus expansion (Magnus, 1954):

$$\begin{array}{}\text{(71)}& {\mathbf{L}}_{\mathrm{eff}}=\underset{\mathrm{0}}{\overset{t}{\int }}\mathbf{L}\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{t}{\int }}\left[\mathbf{L}\left(t_{\mathrm{1}}\right),\mathbf{L}\left(t_{\mathrm{2}}\right)\right]\mathrm{d}{t}_{\mathrm{1}}\mathrm{d}{t}_{\mathrm{2}}+O\left({\mathbf{L}}^{\mathrm{3}}\right).\end{array}$$

Although the convergence radius of this series is rather small (Maricq, 1982), it can nevertheless be used to check whether the effective Liouvillian can be replaced with the time-averaged Liouvillian (first term on the right-hand side of Eq. 71). If the Liouvillian commutes at all times with itself, all higher-order terms vanish and only the zeroth-order term survives. In this case, the effective Liouvillian is equal to the time-averaged Liouvillian. This happens in addition to the case of a constant Liouvillian (see above) if the Liouvillian matrix can be written as a product of a scalar function λ(t) with a constant matrix:

$$\begin{array}{}\text{(72)}& \mathbf{L}\left(t\right)=\mathit{\lambda }\left(t\right)\mathbf{A}\to {\mathbf{L}}_{\mathrm{eff}}=\mathbf{A}\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}.\end{array}$$

In this case, the propagator matrix is exactly the matrix exponential

$$\begin{array}{}\text{(73)}& \mathbf{U}\left(t\right)=\exp \left[-i\mathbf{A}\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}\right].\end{array}$$

This concerns all 2-D cases. Equation (31) becomes

$$\begin{array}{}\text{(74)}& \begin{array}{l}{\mathbf{L}}_{\mathrm{2}\mathrm{D}}=\mathit{\lambda }\left(t\right)\left(\begin{array}{cc}\mathrm{0}& -i\\ i& \mathrm{0}\end{array}\right),\\ \begin{array}{rl}{\mathbf{U}}_{\mathrm{2}\mathrm{D}}& =\exp \left(-i\underset{\mathrm{0}}{\overset{t}{\int }}{\mathbf{L}}_{\mathrm{2}\mathrm{D}}\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}\right)\\ & =\left(\begin{array}{cc}\cos {\int }_{\mathrm{0}}^t\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}& -\sin {\int }_{\mathrm{0}}^t\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}\\ \sin {\int }_{\mathrm{0}}^t\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}& \cos {\int }_{\mathrm{0}}^t\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}\end{array}\right).\end{array}\end{array}\end{array}$$

The transformation to propagation formulae gives

$$\begin{array}{}\text{(75)}& \begin{array}{l}\widehat{A}\stackrel{\widehat{H}t}{⟶}\widehat{A}\cos \underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}+\widehat{B}\sin \underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}},\\ \widehat{B}\stackrel{\widehat{H}t}{⟶}\widehat{B}\cos \underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}-\widehat{A}\sin \underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\lambda }\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}.\end{array}\end{array}$$

Equation (75) can be regarded as a rigorous extension of the POF to time-dependent systems, such as those caused by thermal motion or sample spinning. For higher-dimensional cases, Eq. (72) can only be satisfied if the two parameters a and b are equal. While for the 4-D, 5-D, and 6-D cases this would mean very special situations, the 3-D case includes with a=b the important case of cross-polarization. Then Eq. (38) changes to

$$\begin{array}{}\text{(76)}& {\mathbf{L}}_{\mathrm{3}\mathrm{D}}=a\left(t\right)\left(\begin{array}{ccc}\mathrm{0}& -i& \mathrm{0}\\ i& \mathrm{0}& -i\\ \mathrm{0}& i& \mathrm{0}\end{array}\right)\end{array}$$

if b is replaced with a. In the subsequent propagation formulae, we have to substitute qt with ${\int }_{\mathrm{0}}^{t}q\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}$, with $q\left(t\right)=a\left(t\right)\sqrt{\mathrm{2}}$. This means that we obtain an exact solution if we replace the arguments of the trigonometric functions with integrals in the following examples from above:

• For Eqs. (32), (35), and (52), replace D_IIt with ${\int }_{\mathrm{0}}^t{D}_{\mathrm{II}}\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}$.

• For Eqs. (33), (36), and (53), replace D_ISt with ${\int }_{\mathrm{0}}^t{D}_{\mathrm{IS}}\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}$.

• For Eq. (34), replace ω_Qt with ${\int }_{\mathrm{0}}^t{\mathit{\omega }}_Q\left(t_{\mathrm{1}}\right)\mathrm{d}{t}_{\mathrm{1}}$.

However, there is no general recipe for all the other cases. In some of the cases, the Shirley method using the Floquet theorem, which is applied, e.g., for numerical calculations of spin systems under MAS, may also be successful in obtaining analytical expressions. This will be the subject of a forthcoming paper.

6 Conclusions

Repeated application of the commutator of the Hamiltonian with the initial density operator gives a system of operator equations, the coefficient matrix of which can be used to establish a propagation rule for the spin system. This has been demonstrated in this paper with some examples. A more detailed analysis shows that the commutator relations define subspaces which are both Liouvillian-invariant and superpropagator-invariant. Therefore, the density operator propagates in such a subspace without leaving it. If its dimension is small enough, analytical expressions for the propagation law can be obtained.

The relevant subspace for a given problem is determined by the Hamiltonian and by the initial state. If the operator characterizing the initial state changes, then the new subspace is the same as the previous one if and only if the new initial state operator is an element of the previous subspace. Otherwise, the two subspaces have no intersection.

The set of problems can be divided into classes with respect to the dimension of the subspaces. Problems of the 2-D class can be treated easily by propagation formulae similar to those of the well-known product–operator formalism. The propagation formulae for the 3-D and 4-D classes are given in this paper, and an example is given for the 5-D and 6-D cases. If necessary, the method introduced and explained here can also be applied to cases with higher dimensions. The application examples demonstrate the same mathematical structure as some physically different problems.

In addition, this treatment can be applied to pulse sequences in a manner similar to the POF. In some cases, an algebraic language program may be helpful. Even numerical computations can use this framework by starting the numerical computations on the basis of existing analytical relations. In the N-dimensional wave-function space, the Liouville–von Neumann equation corresponds to a system of N² differential equations, which is significantly larger than that obtained by the dimension reduction with the method presented here. It may be advantageous to first apply this method on an analytical basis before starting the numerical implementation.

If the strength of any of the considered interactions depends on time, the first two steps of this method can be applied likewise. The third step (matrix exponentialization), however, has to be modified because there is no general recipe for solving a system of differential equations with time-dependent coefficients. In some cases, the time-averaged Liouvillian is suitable for getting an exact solution to the problem. Even for numerical calculations it can be helpful to start with systems of differential equations containing a reduced number of equations.

Appendix A: Norms of spin operators

Consider a system of N spins $\mathrm{1}/\mathrm{2}$. As proven in the Supplement (Sect. S1.3), the norm of a product of the Cartesian components of n spins from this N-spin system (one operator for each spin) amounts to ${\mathrm{2}}^{(N/\mathrm{2})-n}$. Furthermore, for the sum of two orthogonal operators $\widehat{A}$ and $\widehat{B}$, $\parallel \widehat{A}+\widehat{B}\parallel =\sqrt{\parallel \widehat{A}{\parallel }^{\mathrm{2}}+\parallel \widehat{B}{\parallel }^{\mathrm{2}}}$. From this, we can derive the following rules (α, β, γ, and $\mathit{\delta }\in \mathit{\{}x,y,z\mathit{\}}$):

• For single spins $I=\frac{\mathrm{1}}{\mathrm{2}}$:
  $\parallel {\widehat{I}}_{\mathit{\alpha }}\parallel =\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}$.

• For spin pairs $I=\frac{\mathrm{1}}{\mathrm{2}}$, $S=\frac{\mathrm{1}}{\mathrm{2}}$:
  $\parallel {\widehat{I}}_{\mathit{\alpha }}\parallel =\parallel {\widehat{S}}_{\mathit{\alpha }}\parallel =\mathrm{1}$, $\parallel {\widehat{I}}_{\mathit{\alpha }}{\widehat{S}}_{\mathit{\beta }}\parallel =\frac{\mathrm{1}}{\mathrm{2}},\parallel {\widehat{I}}_{\mathit{\alpha }}{\widehat{S}}_{\mathit{\beta }}\pm {\widehat{I}}_{\mathit{\gamma }}{\widehat{S}}_{\mathit{\delta }}\parallel =\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}$.

• For spin triples $I=\frac{\mathrm{1}}{\mathrm{2}}$, $J=\frac{\mathrm{1}}{\mathrm{2}}$, $S=\frac{\mathrm{1}}{\mathrm{2}}$:
  $\parallel {\widehat{I}}_{\mathit{\alpha }}\parallel =\parallel {\widehat{J}}_{\mathit{\alpha }}\parallel =\parallel {\widehat{S}}_{\mathit{\alpha }}\parallel =\sqrt{\mathrm{2}}$, $\parallel {\widehat{I}}_{\mathit{\alpha }}{\widehat{S}}_{\mathit{\beta }}\parallel =\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}$,
  $\parallel {\widehat{I}}_{\mathit{\alpha }}{\widehat{S}}_{\mathit{\beta }}\pm {\widehat{I}}_{\mathit{\gamma }}{\widehat{S}}_{\mathit{\delta }}\parallel =\mathrm{1}$, $\parallel {\widehat{I}}_{\mathit{\alpha }}{\widehat{J}}_{\mathit{\beta }}{\widehat{S}}_{\mathit{\gamma }}\parallel =\frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathrm{2}}}$.

Similarly, for a single spin, I=1 can be obtained:

$$\parallel {\widehat{I}}_{\mathit{\alpha }}\parallel =\sqrt{\mathrm{2}},\parallel {\widehat{I}}_{\mathit{\alpha }}{\widehat{I}}_{\mathit{\beta }}+{\widehat{I}}_{\mathit{\beta }}{\widehat{I}}_{\mathit{\alpha }}\parallel =\sqrt{\mathrm{2}}.$$

Appendix B: Hamiltonians used in the main part and in the Supplement

• Interaction of the I and S spins with the rf field, i.e., strengths ω_1I and ω_1S, respectively, assumed to be constant and parallel to the x axis in the rotating frame:

  $$\begin{array}{}\text{(B1)}& {\widehat{H}}_{Ix}=-{\mathit{\omega }}_{\mathrm{1}I}{\widehat{I}}_x;{\widehat{H}}_{Sx}=-{\mathit{\omega }}_{\mathrm{1}S}{\widehat{S}}_x.\end{array}$$

  Unless indicated otherwise, the irradiation occurs at the Larmor frequencies of the respective spins. The rf phase is always such that B₁ is parallel to the x axis of the rotating frame, but the calculations can easily be modified to include other directions of the rf field.

• Homonuclear dipolar interaction between spin I₁ and spin I₂ (secular part):

  $$\begin{array}{}\text{(B2)}& {\widehat{H}}_{\mathrm{II}}=D_{\mathrm{II}}\left({\widehat{\mathit{I}}}_{\mathrm{1}}{\widehat{\mathit{I}}}_{\mathrm{2}}-\mathrm{3}{\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}z}\right),\end{array}$$

  where $D_{\mathrm{II}}=\left({\mathit{\mu }}_{\mathrm{0}}{\mathit{\gamma }}_I^{\mathrm{2}}\mathrm{\hslash }/\left(\mathrm{4}\mathit{\pi }{r}_{\mathrm{II}}^{\mathrm{3}}\right)\right)\left(\mathrm{3}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{II}}-\mathrm{1}\right)$, r_II is the length of the vector connecting both spins, θ_II is the angle of this vector with the external magnetic field, γ_I is the gyromagnetic ratio of the I spins, and μ₀ is the permeability of the vacuum.

• Heteronuclear dipolar interaction between spin I and spin S (secular part):

  $$\begin{array}{}\text{(B3)}& {\widehat{H}}_{\mathrm{IS}}=-D_{\mathrm{IS}}\cdot \mathrm{2}{\widehat{I}}_z{\widehat{S}}_z,\end{array}$$

  where $D_{\mathrm{IS}}=\frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}\frac{\mathrm{\hslash }{\mathit{\gamma }}_I{\mathit{\gamma }}_S}{r_{\mathrm{IS}}^{\mathrm{3}}}\left(\mathrm{3}{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{IS}}-\mathrm{1}\right)$, r_IS is the length of the vector connecting both spins, θ_IS is the angle of this vector with the external magnetic field, and γ_S is the gyromagnetic ratio of the S spins. The case of J coupling between nonequal spins is mathematically equivalent to this; the corresponding formulae can be obtained by replacing D_IS with −πJ.

• Matched Hartmann–Hahn cross-polarization between spin I and spin S in the doubly rotating frame (Hartmann and Hahn, 1962),

  $$\begin{array}{}\text{(B4)}& {\widehat{H}}_{\mathrm{HH}}=-D_{\mathrm{IS}}\cdot \left({\widehat{I}}_x{\widehat{S}}_x+{\widehat{I}}_y{\widehat{S}}_y\right),\end{array}$$

  assuming D_IS≪ω_1I and ω_1S so that rapidly oscillating terms can be neglected.

• First-order interaction of the nuclear quadrupole moment with the electric field gradient (secular part):

  $$\begin{array}{}\text{(B5)}& {\widehat{H}}_Q={\mathit{\omega }}_Q\cdot \left({\widehat{I}}_z^{\mathrm{2}}-{\displaystyle \frac{I(I+\mathrm{1})}{\mathrm{3}}}\right),\end{array}$$

  where ω_Q is the quadrupolar frequency.

• Resonance offset Δω in the rotating frame:

  $$\begin{array}{}\text{(B6)}& {\widehat{H}}_{\mathrm{\Delta }}=-\mathrm{\Delta }\mathit{\omega }{\widehat{I}}_z.\end{array}$$

The coupling frequencies are assumed to be unique, i.e., with no powder-averaging unless otherwise stated.