# A continuous approach to Floquet theory
for pulse-sequence optimization in solid-state NMR

Matías Chávez  Physical Chemistry, ETH Zürich, Vladimir-Prelog-Weg 2, 8093
Zürich, Switzerland  Matthias Ernst<EMAIL_ADDRESS>Physical Chemistry, ETH
Zürich, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland

###### Abstract

We present a framework that uses a continuous frequency space to describe and
design solid-state NMR experiments. The approach is similar to the well
established Floquet treatment for NMR, but is not restricted to periodic
Hamiltonians and allows the design of experiments in a reverse fashion. The
framework is based on perturbation theory on a continuous Fourier space, which
leads to effective, i.e., time-independent, Hamiltonians. It allows the back
calculation of the pulse scheme from the desired effective Hamiltonian as a
function of spin-system parameters. We show as an example how to back
calculate the rf irradiation in the MIRROR experiment from the desired
chemical-shift offset behaviour of the sequence.

††preprint: AIP/123-QED

##  I Introduction

Time-dependent Hamiltonians are very common in (solid-state) NMR due to sample
rotation, e.g., magic angle spinning (MAS) Andrew, Bradbury, and Eades (1958,
1959); Lowe (1959) or even rotations about multiple axis as implemented in
double rotation (DOR) Samoson, Lippmaa, and Pines (1988) or dynamic-angle
spinning (DAS) Mueller  _et al._ (1990), and the application of pulse
sequences that can be described by interaction-frame transformations.Haeberlen
(1976); Mehring (1983) Calculating the time evolution of the density operator
under a time-dependent Hamiltonian is usually done numerically using time-
slicing methods where we assume that the Hamiltonian is time constant for
small enough time periods, $\tau$. However, in many cases an analytical
solution to the Liouville-von Neumann equation is advantageous especially if
predictions of the response of the sequence on spin-system parameters are of
interest. Since the late 1960s, various methods have been developed to
approximate such time-dependent Hamiltonians using time-independent
representations. This can be achieved using average Hamiltonian theory (AHT),
Haeberlen and Waugh (1968); Haeberlen (1976); Maricq (1982); Ernst,
Bodenhausen, and Wokaun (1990) Floquet theory Floquet (1883); Shirley (1965);
Scholz, van Beek, and Ernst (2010); Leskes, Madhu, and Vega (2010); Ivanov
_et al._ (2021) or other less-established methods like the Fer expansion Madhu
and Kurur (2006); Takegoshi  _et al._ (2015), the Floquet-Magnus expansion
Mananga and Charpentier (2016, 2018), or the path-sum method. Giscard and
Bonhomme (2020)
Average Hamiltonian theory Haeberlen and Waugh (1968); Haeberlen (1976);
Maricq (1982); Ernst, Bodenhausen, and Wokaun (1990) can be used for periodic
Hamiltonians and is based on the Magnus expansion. Blanes  _et al._ (2009) It
generates a series of time-independent average Hamiltonians that approximate
the time evolution of the density operator over the basic time period with
increasing accuracy. Limitations of AHT are the stroboscopic observations at
integer multiples of the basic time period, e.g., the difficulty to describe
spinning side bands in MAS spectra, and the difficulty to describe
Hamiltonians with multiple incommensurate frequencies that appear for example
in homonuclear Mote, Agarwal, and Madhu (2016) or heteronuclear Ernst (2003)
decoupling experiments. Despite these limitations, AHT is one of the most used
techniques in NMR to develop and optimize pulse sequences.
In Floquet theory Floquet (1883); Shirley (1965); Scholz, van Beek, and Ernst
(2010); Leskes, Madhu, and Vega (2010); Ivanov  _et al._ (2021) the finite-
dimensional time-dependent spin-Hilbert space Hamiltonian is replaces by an
infinite-dimensional but time-independent Floquet Hamiltonian. Both
representations are equivalent and describe the identical time evolution of
the density operator. Floquet theory can be used to describe periodic time-
dependent Hamiltonians without the requirement of stroboscopic sampling, i.e.,
side bands under MAS are predicted correctly Ivanov  _et al._ (2021) and also
multiple incommensurate time dependencies can be included using multi-mode
Floquet theory. Scholz, van Beek, and Ernst (2010); Leskes, Madhu, and Vega
(2010); Ivanov  _et al._ (2021) Using the time-independent Floquet Hamiltonian
makes no approximations but it is not always easy to get physical insights
from the infinite-dimensional matrices. Therefore, operator-based Floquet
theory Augustine (1995); Boender (1996); Boender, Vega, and de Groot (1998);
Ernst, Samoson, and Meier (2005); Ramachandran and Griffin (2005); Scholz,
Meier, and Ernst (2007), that is based on the van Vleck perturbation treatment
Van Vleck (1929); Primas (1961, 1963) in the Floquet space was developed. Such
a treatment allows the analytical block diagonalization of the Floquet
Hamiltonian and a subsequent projection back into the spin-Hilbert space,
generating a series of effective Hamiltonians that describe the time evolution
of the density operator. Floquet theory and especially the operator-based
implementation has become an important tool for the understanding of magnetic-
resonance experiments and also for the design of new experiments.Ivanov  _et
al._ (2021)
In this publication, we present a modified version of Floquet theory that is
based on a continuous frequency space and not a discrete Fourier series. We
show that the modified approach allows the description of sequences of limited
length and also allows, in some cases, a simpler back calculation of pulse-
sequence parameters.

##  II Theory

###  II.1 Standard Floquet theory

The formalism presented in this article is similar to the formulation of
Floquet theory used in magnetic resonance. Therefore we first want to give a
short review of the Floquet treatment and subsequently present the new
formalism. A more detailed derivation can be found in the literature. Shirley
(1965); Leskes, Madhu, and Vega (2010); Ivanov  _et al._ (2021) The time
dependence of the density operator is given by the Liouville-von Neumann
equation

$\displaystyle\dot{\sigma}(t)=-i\left[{\mathcal{H}}(t),\sigma(t)\right]$ (1)

with the solution

$\displaystyle\sigma(t)=\textrm{U}(t)\sigma(0)\textrm{U}^{-1}(t)$ (2)

where the propagator is formally given by

$\displaystyle\textrm{U}(t)={\mathcal{T}}\exp\left(-i\int\limits_{0}^{t}{\mathcal{H}}(\tau)d\tau\right)$
(3)

and $\mathcal{T}$ is the Dyson time-ordering operator. Dyson (1949) An
equivalent formulation and a good starting point for the derivation of Floquet
theory is the differential equation in the propagator

$\displaystyle\dot{\textrm{U}}(t)=-i\,{\mathcal{H}}(t)\textrm{U}(t)$ (4)

Floquet theory for magnetic resonance is based on the assumption that the
Hamiltonian is periodic and, therefore, the Schrödinger equation is a periodic
differential equation. We can incorporate the periodicity of the Hamiltonian
explicitly by expanding it in a Fourier series as

$\displaystyle{\mathcal{H}}(t)=\sum_{k}{\mathcal{H}}^{(k)}e^{ik\omega t}$ (5)

The form of the solution of the periodic differential Eq. (4) is given by the
Floquet theorem Floquet (1883) as

$\displaystyle\textrm{U}(t)=u(t)e^{i\Lambda t}u^{-1}(0)$ (6)

where $\Lambda$ is a diagonal real matrix because ${\mathcal{H}}(t)$ is a
Hermitian operator. The matrix $u(t)$ has the same periodicity as the
Hamiltonian and can also be expressed as a Fourier series

$\displaystyle u(t)=\sum_{n}u_{n}e^{in\omega t}$ (7)

Notice that now the time dependence of the operators is only in the exponent
and, therefore, inserting Eqs. (5) - (7) into Eq. (4) will lead to the
algebraic equation

$\displaystyle\sum_{k}({\mathcal{H}}^{(n-k)}+\omega\delta_{nk})u_{k}=\Lambda
u_{n}$ (8)

This equation is equivalent to a time-independent Schrödinger equation for the
Fourier coefficients $u_{n}$, where the Hamiltonian

$\displaystyle{\mathcal{H}}_{F}=\sum_{k}{\mathcal{H}}^{(n-k)}+\omega\delta_{nk}$
(9)

is called the Floquet Hamiltonian. It is convenient to represent the Floquet
Hamiltonian and $u_{n}$ in a product space of the spin Hilbert space and the
Fourier space, where each basis state corresponds to a Fourier harmonic.
Typically the basis states are written as
$\ket{n,\mu}=\ket{n}\otimes\ket{\mu}$, where $\ket{\mu}$ denote the basis
states of the spin Hilbert space, and $\ket{n}$ are the basis states of the
Fourier space. Throughout this article Greek letters are used for the spin
basis and Latin letters for the Fourier harmonics. Notice, that there are
infinitely many Fourier harmonics ($n\in\mathbb{Z}$), hence the matrix
representation of the operators has infinite dimensions as well. The explicit
representation of ${\mathcal{H}}_{F}$ Scholz, van Beek, and Ernst (2010);
Leskes, Madhu, and Vega (2010); Ivanov  _et al._ (2021) is given by:

$\displaystyle{\mathcal{H}}_{F}=\sum_{n}F_{n}\otimes{\mathcal{H}}^{(n)}+\omega
F_{z}\otimes\mathbf{1}$ (10)

Here, the $F$ operators act on the Fourier space and are defined by
$F_{z}\ket{n}=n\ket{n}$ and $F_{n}\ket{m}=\ket{m+n}$. The solution of Eqs. (8)
or (10) leads to the formal solution of Eq. (4) expressed in the spin Hilbert
space:

$\displaystyle\textrm{U}(t)=\sum_{n}\bra{n}\exp\left(i{\mathcal{H}}_{F}t\right)\ket{0}e^{in\omega
t}$ (11)

In conclusion, we solved the Liouville-von Neumann equation (Eq. (1)) by
converting it to an algebraic equation using Fourier series expansions and the
Floquet theorem assuming a periodic time-dependent Hamiltonian. As a basis we
chose the product of the spin Hilbert space and Fourier space basis, which
uses Fourier harmonics. Typically van Vleck perturbation theory is applied on
the Floquet Hamiltonian to obtain effective Hamiltonians. The procedure can be
found in the literature Ramesh and Krishnan (2001); Vinogradov, Madhu, and
Vega (2001); Ernst, Samoson, and Meier (2005); Scholz, Meier, and Ernst
(2007); Leskes, Madhu, and Vega (2010); Ivanov  _et al._ (2021) and leads to
the first and second-order Hamiltonian

$\displaystyle\bar{{\mathcal{H}}}^{(1)}={\mathcal{H}}^{(0)}$ (12)

and

$\displaystyle\bar{{\mathcal{H}}}^{(2)}=\frac{1}{2}\sum\limits_{n\neq
0}\frac{[{\mathcal{H}}^{(n)},{\mathcal{H}}^{(-n)}]}{n\omega}$ (13)

Extending operator-based Floquet theory to multiple frequencies is
straightforward. The details can be found in several reviews. Scholz, van
Beek, and Ernst (2010); Leskes, Madhu, and Vega (2010); Ivanov  _et al._
(2021)

###  II.2 Frequency-domain formulation

Let us consider an arbitrary time-dependent Hamiltonian that is not
necessarily periodic in time. In this case we cannot use the Floquet theorem
nor a Fourier series to solve Eq. (1), since both require periodicity.
However, we can use similar methods such as the Fourier transform to obtain an
algebraic equation from Eq. (4). Instead of a Fourier series expansion, we
apply the Fourier transformations to each matrix element of the Hamiltonian
${\mathcal{H}}(t)$ and the propagator $\textrm{U}(t)$. This will lead to a
similar description, where the Fourier transformed operators take the role of
the Fourier coefficients in the Floquet approach. In addition, Fourier
transformation leads to a continuous spectrum of frequencies instead of the
Fourier harmonics.
For the derivation of the generalized framework we again start from the
differential equation Eq. (4). As a first step we use the Fourier
transformation (element wise) to define a frequency-domain Hamiltonian and
propagator as

$\displaystyle{\mathcal{H}}(t)=\mathcal{F}^{-1}\\{{\widehat{{\mathcal{H}}}(\Omega)}\\}=\tau\int\limits_{-\infty}^{\infty}\widehat{{\mathcal{H}}}(\Omega)e^{i\Omega
t}\textrm{d}\Omega$ (14)
$\displaystyle\widehat{{\mathcal{H}}}(\Omega)=\mathcal{F}\\{{{\mathcal{H}}(t)}\\}=\frac{1}{2\,\pi\,\tau}\int\limits_{0}^{\tau}{\mathcal{H}}(t)e^{-i\Omega
t}\textrm{d}t$ (15)
$\displaystyle\textrm{U}(t)=\mathcal{F}^{-1}\\{{\widehat{\textrm{U}}(\Omega)}\\}=\tau\int\limits_{-\infty}^{\infty}\widehat{\textrm{U}}(\Omega)e^{i\Omega
t}\textrm{d}\Omega$ (16)
$\displaystyle\widehat{\textrm{U}}(\Omega)=\mathcal{F}\\{\textrm{U}(t)\\}=\frac{1}{2\,\pi\,\tau}\int\limits_{0}^{\tau}\textrm{U}(t)e^{-i\Omega
t}\textrm{d}t$ (17)

The Fourier transformation is normalized such that it leads to a
correspondence with the standard Floquet approach. Thus, the presented theory
is a generalization of standard Floquet theory and reproduces all of its
results. Notice that the normalization also has to ensure
$\mathcal{F}\\{\mathcal{F}^{-1}\\{f(x)\\}\\}==f(x)$. The ’wide’ hats on
$\widehat{\textrm{U}}(\Omega)$ and $\widehat{{\mathcal{H}}}(\Omega)$ indicate
that these are frequency-domain operators. The underlying frequency domain is
continuous, orthonormal, and complete, hence conceptually similar to the
position or momentum basis used, for example, to describe a free particle.
Inserting the expressions in Eq. (4) and using the convolution theorem leads
to

$\displaystyle i\,\mathcal{F}^{-1}\\{\Omega\,\widehat{\textrm{U}}(\Omega)\\}$
$\displaystyle=-i\,\mathcal{F}^{-1}\\{\widehat{{\mathcal{H}}}(\Omega)\\}\,\mathcal{F}^{-1}\\{\widehat{\textrm{U}}(\Omega)\\}$
$\displaystyle=-i\,\mathcal{F}^{-1}\\{\widehat{{\mathcal{H}}}(\Omega)*\widehat{\textrm{U}}(\Omega)\\}$
(18)

where $*$ symbolizes the convolution operation in frequency space. Applying
Fourier transformation to Eq. (18) results in the time-independent equation

$\displaystyle\Omega\,\widehat{\textrm{U}}(\Omega)=-\widehat{{\mathcal{H}}}(\Omega)*\widehat{\textrm{U}}(\Omega)$
(19)

Since convolution with the delta function is an identity operation we can
write Eq. (19) as

$\displaystyle\Omega\,\delta(\Omega)*\widehat{\textrm{U}}(\Omega)=-\widehat{{\mathcal{H}}}(\Omega)*\widehat{\textrm{U}}(\Omega)$
(20)

and express it as a single integral

$\displaystyle\int\limits_{-\infty}^{\infty}\big{[}\Omega\,\delta(\Omega-\Omega^{\prime})+\widehat{{\mathcal{H}}}(\Omega-\Omega^{\prime})\big{]}\widehat{\textrm{U}}(\Omega^{\prime})\,\textrm{d}\Omega^{\prime}=0$
(21)

Since this has to be fulfilled for every possible
$\widehat{\textrm{U}}(\Omega^{\prime})$ the kernel has to vanish

$\displaystyle\big{[}\Omega\,\delta(\Omega-\Omega^{\prime})+\widehat{{\mathcal{H}}}(\Omega-\Omega^{\prime})\big{]}\widehat{\textrm{U}}(\Omega^{\prime})=0$
(22)

Notice, that this is an eigenvalue equation of
$\widehat{\textrm{U}}(\Omega^{\prime})$ with the eigenoperator

$\displaystyle{\mathcal{H}}_{F}=\Omega\,\delta(\Omega-\Omega^{\prime})+\widehat{{\mathcal{H}}}(\Omega-\Omega^{\prime})$
(23)

Similar to the standard Floquet treatment, we can express this Hamiltonian in
the orthonormal product basis $\ket{\Omega,\mu}=\ket{\Omega}\otimes\ket{\mu}$.
As before the states $\ket{\mu}$ denote the basis states of the spin Hilbert
space, but $\ket{\Omega}$ are the basis states of the continuous Fourier
space. Notice that the operators
$\widehat{{\mathcal{H}}}(\Omega-\Omega^{\prime})$ and $\Omega$ in
${\mathcal{H}}_{F}$ live in different subspaces, i.e the spin and the
frequency space, respectively. Hence the Hamiltonian in the product basis
takes the form

$\displaystyle\bra{\Omega,\mu}{\mathcal{H}}_{F}\ket{\Omega^{\prime},\nu}=\Omega\,\delta(\Omega-\Omega^{\prime})\,\otimes\delta_{\mu\nu}$
$\displaystyle+\mathbf{1}\otimes\widehat{{\mathcal{H}}}_{\mu\nu}(\Omega-\Omega^{\prime})$
(24)

where
$\widehat{{\mathcal{H}}}_{\mu\nu}(\Omega)=\bra{\mu}\widehat{{\mathcal{H}}}(\Omega)\ket{\nu}$
are the matrix elements of the spin-system Hamiltonian at frequency $\Omega$.

In full analogy to the Floquet description (Eq. (10)), ${\mathcal{H}}_{F}$ can
be formulated in a basis-free operator form as

$\displaystyle{\mathcal{H}}_{F}=\int\widehat{D}(\Omega)\otimes\widehat{{\mathcal{H}}}(\Omega)\,\textrm{d}\Omega+\widehat{\Omega}\otimes\mathbf{1}$
(25)

with the frequency operator

$\displaystyle\widehat{\Omega}\ket{\Omega}=\Omega\ket{\Omega}$
$\displaystyle\iff\widehat{\Omega}=\int\Omega^{\prime}\ket{\Omega^{\prime}}\bra{\Omega^{\prime}}\,\textrm{d}\Omega^{\prime}$
(26)

and the frequency-shift operator

$\displaystyle\widehat{D}(\Omega_{0})\ket{\Omega}=\ket{\Omega+\Omega_{0}}$
(27)
$\displaystyle\widehat{D}(\Omega_{0})=\int_{\Omega^{\prime}}\ket{\Omega^{\prime}+\Omega_{0}}\bra{\Omega^{\prime}}\textrm{d}\Omega^{\prime}.$
(28)

The equivalence of Eq. (25) and Eq. (23) is shown in the supplementary
information (SI). Notice that we introduced the frequency operator
$\widehat{\Omega}$, which is similar to a position operator, but acting on the
frequency space. It takes the place of the number operator $F_{z}$ in standard
Floquet theory. The frequency-shift operator $\widehat{D}(\Omega)$ is
conceptually similar to the usual translation operator in the position space.
It takes the place of the generalized ladder operators $F_{n}$ in standard
Floquet theory. Both of this operators act exclusively on the frequency
domain.
As in the standard Floquet treatment used in solid state NMR we utilize a
frequency basis together with the spin basis to represent the Hamiltonian. In
contrast to the standard Floquet theory we use a continuous frequency basis,
which enables the description of an arbitrary time-modulated Hamiltonian.
Nonetheless, we obtained a similar description, employing a continuous Fourier
basis and Fourier transformations instead of the discrete Fourier basis and
Fourier series as in the Floquet approach. Since we do not have a discrete
Fourier basis, we do not obtain a representation of the Hamiltonian as a
matrix of constant coefficients, but as a frequency dependent matrix function.
However, the commutation relations are similar to those of the Floquet
approach and allow for a similar treatment for example using van Vleck
perturbation theory. Van Vleck (1929); Primas (1961, 1963)

####  II.2.1 Generalization to multiple time-dependent modulations

In principle, with the continuous frequency basis, there is no need to
describe multiple time-dependent modulations using different frequency
dimensions as is required in standard Floquet theory. All time-dependent
modulations could be lumped into a single dimension. However, to understand
resonance conditions between different processes, e.g., magic-angle spinning
and radio-frequency irradiation of the spins, it is advantageous to separate
this processes in different dimensions represented by multiple frequency
bases. Let us generalize the approach presented in the previous section to $n$
time-dependent modulations of the Hamiltonian. As before we can describe the
Hamiltonian in a frequency domain by using Fourier transformation, resulting
in the expression

$\displaystyle{\mathcal{H}}(t)=\tau^{n}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}...\int\limits_{-\infty}^{\infty}\widehat{{\mathcal{H}}}(\Omega_{1},\Omega_{2},\cdots,\Omega_{n})$
$\displaystyle e^{i\Omega_{1}t}e^{i\Omega_{2}t}\cdots
e^{i\Omega_{n}t}\,\textrm{d}\Omega_{1}\,\textrm{d}\Omega_{2}\cdots\textrm{d}\Omega_{n}$
(29)

Now the underlying basis is the direct product basis of $n$ frequency bases
and one spin basis $\\{\ket{\Omega_{1},\Omega_{2},\cdots,\Omega_{n},\mu}\\}$.
This approach is similar to multimodal Floquet theory, but is not limited to
periodic modulations of the Hamiltonian. From Eq. (29) we directly obtain the
resonance condition

$\displaystyle\Omega_{1}^{(0)}+\Omega_{2}^{(0)}+\cdots+\Omega_{n}^{(0)}=0$
(30)

Following the same procedure as in the single mode case we obtain a time-
independent Hamiltonian of the form

$\displaystyle{\mathcal{H}}_{F}$
$\displaystyle=\int\int\cdots\int\widehat{D}_{1}(\Omega_{1})\otimes\widehat{D}_{2}(\Omega_{2})\otimes\cdots\otimes\widehat{D}_{n}(\Omega_{n})$
$\displaystyle\quad\otimes\widehat{{\mathcal{H}}}(\Omega_{1},\Omega_{2},\cdots,\Omega_{n})\,\,\textrm{d}\Omega_{1}\textrm{d}\Omega_{2}\cdots\textrm{d}\Omega_{n}$
$\displaystyle\quad\oplus\widehat{\Omega}_{1}\oplus\widehat{\Omega}_{2}\oplus\cdots\oplus\widehat{\Omega}_{n}$
(31)

where $\widehat{D}_{j}(\Omega)$ with $j\in\\{1,2,\cdots,n\\}$ is the
corresponding frequency-translation operator and $\widehat{\Omega}_{j}$ the
frequency operator. It is important to emphasise that in practice most likely
two frequency domains will be sufficient, a spatial modulation, for example
due to sample rotation, and a modulation of the spin system due to rf-field
irradiation. In this case the bimodal approach will be sufficient. However,
for rf-field irradiation addressing different spin species, like proton,
nitrogen or electron spins, a higher modal approach might offer advantages. In
principle, one can also combine the frequency-domain approach in one dimension
(e.g., rf irradiation) with the traditional Floquet approach in a different
dimension (e.g., sample rotation) where the Fourier series is a perfect
description of the time-dependent Hamiltonian (vide infra).

###  II.3 Van Vleck perturbation theory and effective Hamiltonian

In the following, we apply van Vleck perturbation theory Van Vleck (1929);
Primas (1961, 1963) to the Floquet Hamiltonian defined in Eq. (25) with the
goal of obtaining effective Hamiltonians of different orders. The detailed
derivation can be found in the SI and follows the treatment in Ref.Ernst,
Samoson, and Meier (2005) As usual we split the Hamiltonian into two parts

$\displaystyle{\mathcal{H}}_{F}={\mathcal{H}}_{F}^{(0)}+\varepsilon\,{\mathcal{H}}_{F}^{(1)}$
(32)

and apply the van Vleck transformation given by

$\displaystyle
e^{S}\,{\mathcal{H}}_{F}\,e^{S^{\dagger}}=\sum_{m=0}^{\infty}\frac{[S,{\mathcal{H}}_{F}]_{m}}{m!}=\sum_{n=1}^{\infty}\varepsilon^{n}\Lambda_{F}^{(n)}$
(33)

where the nested commutator is defined as
$[S,{\mathcal{H}}_{F}^{(0)}]_{m}=[S,[S,{\mathcal{H}}_{F}^{(0)}]_{m-1}]$ with
$[S,{\mathcal{H}}_{F}^{(0)}]_{0}={\mathcal{H}}_{F}^{(0)}$. We choose $S$, such
that $[\Lambda_{F},{\mathcal{H}}_{F}^{(0)}]=0$ and proceed by expanding $S$ as
a series

$\displaystyle S=\sum_{n=1}^{\infty}\varepsilon^{n}S^{(n)}$ (34)

and inserting it into Eq. (33), which leads to

$\displaystyle[S^{(n)},{\mathcal{H}}_{F}^{(0)}]=\Lambda_{F}^{(n)}-\Phi^{(n)}_{F}$
(35)

with

$\displaystyle\sum_{j=1}^{\infty}\varepsilon^{j}\Phi^{(j)}_{F}=$
$\displaystyle\sum_{m=2}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},{\mathcal{H}}_{F}^{(0)}]_{m}}{m!}$
$\displaystyle+\sum_{m=1}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},\varepsilon\,{\mathcal{H}}_{F}^{(1)}]_{m}}{m!}.$
(36)

Adopting the approach from Primas Primas (1961, 1963) we obtain the formal
solution of Eq. (162) as

$\displaystyle S^{(n)}$
$\displaystyle=\Gamma_{F}^{-1}(\Pi(\Phi^{(n)}_{F})-\Phi^{(n)}_{F})$ (37)
$\displaystyle\Lambda_{F}^{(n)}$ $\displaystyle=\Pi(\Phi^{(n)}_{F})$ (38)

where $\Pi(X)$ is the projection operator and $\Gamma_{F}^{-1}$ the inverse
commutation operator defined in the SI.
Let us calculate the first and second-order effective Hamiltonians for the
single-mode case:

$\displaystyle{\mathcal{H}}^{(0)}_{F}$ $\displaystyle=\widehat{\Omega}$ (39)
$\displaystyle{\mathcal{H}}^{(1)}_{F}$
$\displaystyle=\int\widehat{D}(\Omega)\otimes\widehat{{\mathcal{H}}}(\Omega)\,\textrm{d}\Omega=\Phi^{(1)}_{F}$
(40)

For the first-order effective Hamiltonian we evaluate Eq. (159) leading to

$\displaystyle\bar{{\mathcal{H}}}^{(1)}=\bra{\Omega^{\prime}}\Lambda_{F}^{(1)}\ket{\Omega}=\widehat{{\mathcal{H}}}(0)=\frac{1}{2\,\pi\,\tau}\int\limits_{0}^{\tau}{{\mathcal{H}}}(t)\,\textrm{d}t$
(41)

For the second-order Hamiltonian we use Eqs. (153-159) resulting in

$\displaystyle\bar{{\mathcal{H}}}^{(2)}=\frac{1}{2}\,PV\int\frac{[\widehat{{\mathcal{H}}}(\Omega),\widehat{{\mathcal{H}}}(-\Omega)]}{\Omega}\textrm{d}\Omega$
(42)

We use the Cauchy principal value ($PV$) for the regularization of the
integral, which avoids the integration over the singularity using limits. The
main difference of Eq. (42) to Floquet theory is the integral over the
frequency-domain Hamiltonian, instead of a sum of Fourier coefficients of the
Hamiltonian. This disparity stems from the fact that the underlying Fourier
space is continuous in contrast to the discrete space we use to describe
periodic Hamiltonians.
The derivation of the effective Hamiltonians for multiple frequency dimensions
is quite similar and leads to the first and second-order effective Hamiltonian

$\displaystyle\bar{{\mathcal{H}}}^{(1)}$
$\displaystyle=\int\widehat{{\mathcal{H}}}(\Omega^{(0)}_{1},-\Omega^{(0)}_{1})d\Omega^{(0)}_{1}$
(43)

and

$\displaystyle\bar{{\mathcal{H}}}^{(2)}$
$\displaystyle=\frac{1}{2}PV\int\textrm{d}\Omega_{1}\int\textrm{d}\Omega_{1}^{(0)}\int\textrm{d}\Omega_{2}$
$\displaystyle\quad\frac{[\widehat{{\mathcal{H}}}(\Omega_{1},\Omega_{2}),\widehat{{\mathcal{H}}}(\Omega_{1}^{(0)}-\Omega_{1},-\Omega_{1}^{(0)}-\Omega_{2})]}{\Omega_{1}+\Omega_{2}}$
(44)

where we used the resonance condition

$\displaystyle\Omega^{(0)}_{2}=-\Omega^{(0)}_{1}$ (45)

In fact, we can retrieve the effective Hamiltonians obtained from Floquet
theory, considering a periodic Hamiltonian. Periodicity with period
$\tau_{\textrm{m}}$, causes quantization in the frequency domain with the
frequency $\omega_{\textrm{m}}=2\,\pi/\tau_{\textrm{m}}$. Therefore the
frequency-domain Hamiltonian is only non-zero at the harmonics
$n\,\omega_{\textrm{m}}$. Expressing this fact by a Dirac comb we obtain

$\displaystyle{\mathcal{H}}(t)$
$\displaystyle=\int\limits_{-\infty}^{\infty}\widehat{{\mathcal{H}}}(\Omega)\,e^{i\Omega
t}\,\textrm{d}\Omega$
$\displaystyle=\int\limits_{-\infty}^{\infty}\sum_{n=0}^{\infty}\delta(\Omega-n\,\omega_{\textrm{m}})\,\widehat{{\mathcal{H}}}(\Omega)\,e^{i\Omega
t}\,\textrm{d}\Omega$
$\displaystyle=\sum_{n=0}^{\infty}\widehat{{\mathcal{H}}}(n\,\omega_{\textrm{m}})\,e^{in\omega_{\textrm{m}}t}$
$\displaystyle=\sum_{n=0}^{\infty}{\mathcal{H}}^{(n)}\,e^{in\omega_{\textrm{m}}t}$
(46)

where we can identify the Fourier coefficients of the Hamiltonian as
$\widehat{{\mathcal{H}}}(n\,\omega_{\textrm{m}})={\mathcal{H}}^{(n)}$. This is
simply the transition of the Fourier transformation to a Fourier series for a
periodic function. As a consequence imposing a periodic boundary condition on
the Hamiltonian leads to the results from Floquet theory.
As mentioned already above, it is possible to use a mixed approach, describing
periodic modulations with a discrete Fourier space and non-periodic
modulations on a continuous Fourier space. We can use the same approach with
the delta comb as in Eq. (46) to arrive at

$\displaystyle{\mathcal{H}}(t)$
$\displaystyle=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\sum_{n}\delta(\Omega_{1}-n\omega_{\textrm{m}})\widehat{{\mathcal{H}}}(\Omega_{1},\Omega_{2})\,e^{i\Omega_{1}t}e^{i\Omega_{2}t}\,\textrm{d}\Omega_{1}\,\textrm{d}\Omega_{2}$
$\displaystyle=\int\limits_{-\infty}^{\infty}\sum_{n}\widehat{{\mathcal{H}}}(n\omega_{\textrm{m}},\Omega_{2})\,e^{in\omega_{\textrm{m}}t}e^{i\Omega_{2}t}\,\textrm{d}\Omega_{2}$
$\displaystyle=\int\limits_{-\infty}^{\infty}\sum_{n}\widehat{{\mathcal{H}}}^{(n)}(\Omega_{2})\,e^{in\omega_{\textrm{m}}t}e^{i\Omega_{2}t}\,\textrm{d}\Omega_{2}$
(47)

This approach is useful to analyze general solid-state NMR experiments under
MAS, where the sample spinning is described as a periodic modulation, but the
rf-field irradiation can be arbitrary. Such mixed approaches can be used for
any number of modes, where for example the periodic modulation due to the
chemical shift offset is treated as a third mode. The first and second-order
effective Hamiltonians for the mixed approach with two frequencies as in Eq.
(47) have the form

$\displaystyle\bar{{\mathcal{H}}}^{(1)}$
$\displaystyle=\sum_{n_{0}}\sum_{\Omega_{0}}\widehat{{\mathcal{H}}}^{(n_{0})}(\Omega_{0})$
(48)

and

$\displaystyle\bar{{\mathcal{H}}}^{(2)}$
$\displaystyle=\frac{1}{2}\sum_{n}\sum_{n_{0},\Omega_{0}}PV\int\textrm{d}\Omega\quad\frac{[\widehat{{\mathcal{H}}}^{(n)}(\Omega),\widehat{{\mathcal{H}}}^{(n_{0}-n)}(\Omega_{0}-\Omega)]}{\Omega+n\omega_{\textrm{m}}}$
(49)

with the resonance conditions defined by
$\Omega_{0}+n_{0}\omega_{\textrm{m}}=0$. Similar to standard Floquet theory,
the propagation with an effective Hamiltonian is mediated by

$\displaystyle\bar{\textrm{U}}(t)=\exp(-i\,\bar{{\mathcal{H}}}^{(1)}\,t)$ (50)

In contrast to standard Floquet Theory, the effective Hamiltonian is only
valid for one specific duration, since the time dependence is implicitly in
${\mathcal{H}}^{(1)}$. However, if the duration of the irradiation scheme is
changed or the scheme is applied repeatedly, the effective Hamiltonian is
weighted by a function, which depends on the duration or repetitions. In
section II.5.2 we derive this function and discuss its effects.

###  II.4 Calculation of the frequency-domain interaction-frame trajectory

As already mentioned, sample spinning in NMR leads to a periodic modulation of
the Hamiltonian and can be described by a discrete Fourier series with the
Fourier coefficients typically limited to the range -2 to 2. Radio-frequency
irradiation is also often periodic in time, but typically an interaction-frame
transformation is required to ensure the convergence of the effective
Hamiltonian series. We split the Hamiltonian into two parts

$\displaystyle{\mathcal{H}}(t)={\mathcal{H}}_{\textrm{S}}(t)+{\mathcal{H}}_{\textrm{CS}}(t)+{\mathcal{H}}_{\textrm{rf}}(t)={\mathcal{H}}_{\textrm{0}}(t)+{\mathcal{H}}_{1}(t)$
(51)

where ${\mathcal{H}}_{\textrm{0}}$ describes the spin-system Hamiltonian,
${\mathcal{H}}_{\textrm{rf}}(t)$ is the rf-field Hamiltonian and
${\mathcal{H}}_{\textrm{CS}}(t)$ the isotropic chemical-shift Hamiltonian. We
can choose how to divide the total Hamiltonian into ${\mathcal{H}}_{0}$ and
${\mathcal{H}}_{1}$. Either we can set
${\mathcal{H}}_{1}(t)={\mathcal{H}}_{\textrm{rf}}(t)$ and
${\mathcal{H}}_{0}(t)={\mathcal{H}}_{\textrm{S}}(t)+{\mathcal{H}}_{\textrm{CS}}(t)$
resulting in an interaction-frame transformation by the rf irradiation only.
In this case, the interaction-frame transformation is the same for all spins.
Alternatively we can set
${\mathcal{H}}_{1}(t)={\mathcal{H}}_{\textrm{CS}}(t)+{\mathcal{H}}_{\textrm{rf}}(t)$
and ${\mathcal{H}}_{0}(t)={\mathcal{H}}_{\textrm{S}}(t)$ resulting in an
interaction-frame transformation by the time-dependent effective field of each
spin. In this case we have again two options. We can either include only the
isotropic chemical shift into the interaction-frame transformation or we can
include the isotropic and the anisotropic chemically shift. With this choice,
the interaction-frame trajectory of each spin will be distinct if the chemical
shifts are different. The selection of the most convenient interaction frame
will depend on the system and problem at hand. In the following we use an
interaction frame generated by a general rf-field modulation including the
chemical shift offset. In this case the unitary interaction-frame
transformation of the Hamiltonian for each spin is a general complex rotation.
Again, each spin will have its unique unitary transformation and, therefore,
its unique frame, except if they are chemical equivalent. In the usual
rotating frame defined by the Zeeman Hamiltonians, the Hamiltonian used for
the interaction-frame transformation has the form

$\displaystyle{\mathcal{H}}_{1}(t)$
$\displaystyle=\omega_{1}(t)[\cos(\phi(t))\textrm{I}_{x}+\sin(\phi(t))\textrm{I}_{y}]+\omega_{z}(t)\textrm{I}_{z}$
$\displaystyle=\vec{\theta}(t)\vec{\textrm{I}}$ (52)

with

$\displaystyle\vec{\theta}(t)=\begin{pmatrix}\theta_{x}(t)\\\ \theta_{y}(t)\\\
\theta_{z}(t)\\\
\end{pmatrix}=\begin{pmatrix}\omega_{1}(t)\cos\left(\phi(t)\right)\\\
\omega_{1}(t)\sin\left(\phi(t)\right)\\\ \omega_{z}(t)\end{pmatrix}$ (53)
$\displaystyle|\theta(t)|=\sqrt{\omega^{2}_{1}(t)+\omega^{2}_{z}(t)}$ (54)

where $\omega_{1}(t)$ is the amplitude and $\phi(t)$ the phase of the rf-field
and $\omega_{z}(t)$ the chemical shift offset. The direction of
$\vec{\theta}(t)$ specifies the axis of rotation, its length $|\theta(t)|$ is
the angular velocity at time $t$. The interaction-frame transformation is
given by

$\displaystyle\textrm{U}(t)={\cal{T}}\exp\left(-i\,\int\limits_{0}^{t}{\mathcal{H}}_{1}(t^{\prime})\,\textrm{d}t^{\prime}\right)$
(55)

where $\cal{T}$ is the Dyson time-ordering operator. Dyson (1949) The
evolution of the initial spin operators
$\vec{\textrm{I}}=(\textrm{I}_{x},\textrm{I}_{y},\textrm{I}_{z})^{T}$ can be
written with a single rotation matrix with elements $a_{\mu\nu}(t)$ as

$\displaystyle\textrm{I}_{\mu}(t)=\textrm{U}^{\dagger}(t)\,\textrm{I}_{\mu}\,\textrm{U}(t)=\sum\limits_{\nu}a_{\mu\nu}(t)\,\textrm{I}_{\nu}$
(56)

This also can be expressed in the frequency domain using Fourier transform
leading to

$\displaystyle\textrm{I}_{\mu}(t)=\mathcal{F}^{-1}\\{\widehat{\textrm{I}}_{\mu}(\Omega)\\}=\sum\limits_{\nu}a_{\mu\nu}(t)\,\textrm{I}_{\nu}$
(57)
$\displaystyle\widehat{\textrm{I}}_{\mu}(\Omega)=\mathcal{F}\\{\textrm{I}_{\mu}(t)\\}=\sum\limits_{\nu}\widehat{a}_{\mu\nu}(\Omega)\,\textrm{I}_{\nu}$
(58)

where $\widehat{a}_{\mu\nu}(\Omega)=\mathcal{F}\\{a_{\mu\nu}(t)\\}$. In the
following we present an efficient route to calculate $a_{\mu\nu}(t)$. Notice
that the map
$\textrm{I}_{\mu}\mapsto\textrm{U}(t)\,\textrm{I}_{\mu}\,\textrm{U}^{\dagger}(t)$
and
$\textrm{I}_{\mu}\mapsto(-\textrm{U}(t))\,\textrm{I}_{\mu}\,(-\textrm{U}^{\dagger}(t))$
leads to the same rotation of $\mathbb{R}^{3}$, hence two elements of $SU(2)$
are mapped onto one element of $SO(3)$. In more technical terms this means
that there exists a 2:1 surjective homomorphism from $SU(2)$ to $SO(3)$.
However Eq. (55) yields to only one of the two $SU(2)$ elements, since it is
impossible to get $-\textrm{U}(t)$ from a given $\textrm{U}(t)$ with the
Hamiltonian of the form given in Eq. (52). Therefore, the operator defined in
Eq. (55) is always a member of $SU(2)/Z_{2}$ which is isomorphic to $SO(3)$,
i.e. $SU(2)/Z_{2}\cong SO(3)$. This fact will be important for the back-
calculation, when we infer the rf-irradiation from rotation matrices.A second
important point is, that in general, ${\mathcal{H}}_{1}(t)$ does not commute
with itself at different times. For this reason we divide the irradiation into
time intervals, small enough to assume that the Hamiltonian commutes with
itself during the interval. However, the intervals do not necessary have to be
of the same size and can be adjusted to fit the problem at hand. The
propagator of the $j\,$th interval corresponds to $\textrm{U}_{j}$ and has a
complex $2\times 2$ matrix representation given in the SI. Equipped with this
matrix representation we calculate the coefficients of the rotation matrices
$a^{(j)}_{\mu\nu}$ corresponding to $\textrm{U}_{j}$ as

$\displaystyle\textrm{U}_{j}^{\dagger}\,\textrm{I}_{\mu}\,\textrm{U}_{j}=\sum_{\nu}a_{\mu\nu}^{(j)}\textrm{I}_{\nu}$
(59)

Using
$a_{\mu\nu}=2\,\text{Tr}(\textrm{I}_{\mu}\,\textrm{U}\,\textrm{I}_{\nu}\,\textrm{U}^{-1})$,
leads to Cornwell (1984)

$\displaystyle
a_{xx}=\left[\left(\theta^{2}_{y}+\theta_{z}^{2}\right)\cos\left(|\theta|\right)+\theta^{2}_{x}\right]/(4|\theta|)$
$\displaystyle
a_{xy}=\left[-\theta_{x}\theta_{y}\left(\cos\left(|\theta|\right)-1\right)-\theta_{z}|\theta|\sin\left(|\theta|\right)\right]/(4|\theta|)$
(60) $\displaystyle
a_{xz}=\left[\theta_{y}|\theta|\sin\left(|\theta|\right)+\theta_{x}\theta_{z}\left(1-\cos\left(|\theta|\right)\right)\right]/(4|\theta|)$
$\displaystyle
a_{yx}=\left[\theta_{z}|\theta|\sin\left(|\theta|\right)-\theta_{x}\theta_{y}\left(\cos\left(|\theta|\right)-1\right)\right]/(4|\theta|)$
$\displaystyle
a_{yy}=\left[\left(\theta^{2}_{x}+\theta_{z}^{2}\right)\cos\left(|\theta|\right)+\theta^{2}_{y}\right]/(4|\theta|)$
(61) $\displaystyle
a_{yz}=\left[-\theta_{y}\,\theta_{z}\left(\cos\left(|\theta|\right)-1\right)-\theta_{x}|\theta|\sin\left(|\theta|\right)\right]/(4|\theta|)$
$\displaystyle
a_{zx}=\left[-\theta_{y}|\theta|\sin\left(|\theta|\right)-\theta_{z}\theta_{x}\left(\cos\left(|\theta|\right)-1\right)\right]/(4|\theta|)$
$\displaystyle
a_{zy}=\left[\theta_{x}|\theta|\sin\left(|\theta|\right)-\theta_{z}\theta_{y}\left(\cos\left(|\theta|\right)-1\right)\right]/(4|\theta|)$
(62) $\displaystyle
a_{zz}=\left[(\theta^{2}_{x}+\theta^{2}_{y})\cos\left(|\theta|\right)+\theta_{z}^{2}\right]/(4|\theta|)$

where we omitted the indices $j$ and $\theta$ on $a_{\mu\nu}$ for the sake of
simplicity. Equivalent expressions can also be found for the alternate basis
with $\mu,\nu\in\\{+,-,z\\}$. The interaction-frame trajectory at time $t_{j}$
is given by

$\displaystyle{\bf{a}}(t_{j})={\bf{a}}_{j}\,{\bf{a}}_{j-1}\cdots{\bf{a}}_{1}$
(63)

Note, that the rotation matrices ${\bf{a}}_{j}$ describe the rotation during a
single time step while the matrix ${\bf{a}}(t_{j})$ describes the total
rotation up to the time point $t_{j}$. The expressions for the elements
$a_{\mu\nu}(t_{j})$ in dependence on the parameter $\omega_{1}$, $\omega_{z}$,
and $\phi$ can be found in the SI. As mentioned before, the elements
$\widehat{a}_{\mu\nu}{({\Omega})}$ are obtained by Fourier transformation of
the elements $a_{\mu\nu}(t)$

$\displaystyle\widehat{a}_{\mu\nu}(\Omega)=\mathcal{F}\\{a_{\mu\nu}(t)\\}=\frac{1}{2\pi\tau}\int\limits_{0}^{\tau}a_{\mu\nu}(t)e^{i\Omega
t}\,dt$ (64)

The mathematical properties of $a_{\mu\nu}(t)$ and
$\widehat{a}_{\mu\nu}(\Omega)$ can be summarized as follows

$\displaystyle a_{\mu\nu}(t)\in\mathbb{R}$
$\displaystyle\Leftrightarrow\widehat{a}_{\mu\nu}(\Omega)=\widehat{a}^{*}_{\nu\mu}(-\Omega)$
(65) $\displaystyle\sum_{\chi}a_{\mu\chi}(t)a_{\nu\chi}(t)=\delta_{\mu\nu}$
$\displaystyle\Leftrightarrow\sum\limits_{\chi}\widehat{a}_{\mu\chi}(\Omega)*\widehat{a}_{\nu\chi}(\Omega)=\delta_{\mu\nu}\delta(\Omega)$
(66)

###  II.5 Properties of the interaction-frame trajectory

In the following section the properties of the time and frequency-domain
interaction-frame trajectory are explored. We will show that any finite
irradiation can be expressed in term of its periodic version, on one hand
reducing the computational effort, on the other hand isolating the effect of
the finite duration. Furthermore we will show that for any cyclic irradiation
the calculation of the frequency-domain interaction-frame trajectory can be
reduced to the calculation of a single segment, further lessening the
computational efforts significantly. Finally we will discuss the case of
irradiation along a single axis, which always enables the expression of the
frequency-domain interaction-frame trajectory in a closed form.

####  II.5.1 General properties of the interaction-frame trajectory

From the Plancherel theorem we obtain the general property Plancherel and
Leffler (1910)

$\displaystyle\frac{1}{2\pi\tau}\int\limits_{-\infty}^{\infty}|\widehat{a}_{\mu\nu}(\Omega)|^{2}\,d\Omega=\tau\int\limits_{0}^{\tau}|\widehat{a}_{\mu\nu}(t)|^{2}\,dt$
(67)

In contrast to standard Floquet Theory, the duration of the irradiation scheme
is incorporated implicitly in the frequency-domain interaction-frame
trajectory. However any finite $a_{\mu\nu}(t)$ defined on $t\in[0,T]$ can be
made periodic ($t\in\mathbb{R}$) using the modulo function

$\displaystyle\widetilde{a}_{\mu\nu}(t)=a_{\mu\nu}(t\,\textrm{mod}\,T)$ (68)

where we indicate the periodized function using a tilde. For a periodic
interaction-frame trajectory the frequency-domain interaction-frame trajectory
is related to the Fourier coefficients

$\displaystyle\widehat{a}_{\mu\nu}(\Omega)=\sum_{n=-\infty}^{\infty}a_{\mu\nu}^{(n)}\delta(\Omega-n\,{2\pi}/{T})$
(69)

This relation is the bridge between standard Floquet theory, which utilizes
Fourier series and continuous Floquet theory, which uses the Fourier
transform. This correspondence only exists for the periodic case where
standard Floquet theory is valid. Vice versa, a periodic function
$\widetilde{a}_{\mu\nu}(t)$ can be made finite by multiplication with the
appropriate rectangular window function $\Pi(t)$

$\displaystyle
a_{\mu\nu}(t)=\widetilde{a}_{\mu\nu}(t)\cdot\Pi\left(t/T-1/2\right)$ (70)
$\displaystyle\Pi\left(t/T-1/2\right)=\begin{cases}1&0\leq t\leq T\\\ 0&\hskip
3.00003pt\text{else}\end{cases}$ (71)

Eq. (70) and the convolution theorem lead to they frequency-domain
interaction-frame trajectory

$\displaystyle\widehat{a}_{\mu\nu}(\Omega)$
$\displaystyle=\widetilde{a}_{\mu\nu}(\Omega)*\mathcal{F}\left\\{\Pi\left(t/T-1/2\right)\right\\}$
$\displaystyle=\widetilde{a}_{\mu\nu}(\Omega)*\frac{\textrm{sin}(\Omega\,T/2)}{\pi\,\Omega\,T}e^{-i\frac{\Omega}{2}T}$
(72)

Notice, that Eq. (II.5.1) can be efficiently implemented numerically via
discrete Fourier transformation (DFT), since $\widetilde{a}_{\mu\nu}(\Omega)$
is periodic. Furthermore, the slope of $\widehat{a}_{\mu\nu}(\Omega)$ is
dictated by $T$, as $\widetilde{a}_{\mu\nu}(\Omega)$ has sharp frequency
components, see Eq. (69). As a consequence sharp edges of
$\widehat{a}_{\mu\nu}(\Omega)$ cannot be realized by any irradiation with
finite duration. We can go a step further and decompose
$\widetilde{a}_{\mu\nu}(t)$ into a product and subsequently use the
convolution theorem

$\displaystyle\widehat{a}_{\mu\nu}(\Omega)$
$\displaystyle=\mathcal{F}\left\\{\Pi(t/T)\cdot\prod\limits_{i}^{n}\,\widetilde{a}^{(i)}_{\mu\nu}(t)\right\\}$
$\displaystyle=\mathcal{F}\left\\{\Pi(t/T)\right\\}*\widetilde{a}^{(1)}_{\mu\nu}(\Omega)*\widetilde{a}^{(2)}_{\mu\nu}(\Omega)*\ldots*\widetilde{a}^{(n)}_{\mu\nu}(\Omega)$
(73)

This trick will be useful to construct pulse schemes, since we can decompose
the desired $\widehat{a}_{\mu\nu}(\Omega)$ with convolutions and calculate the
necessary $a^{(j)}_{\mu\nu}$ for each part. Finally we obtain the desired
trajectory $a_{\mu\nu}(t)$ by calculating the product
$\prod^{n}_{j}a^{(j)}_{\mu\nu}$. The ability to describe sequences of
different length and the consequences of finite recoupling sequences is one of
the advantages of the new method. Independent of the basis of the spin
operator, $a_{\mu\nu}(t)$ are the elements of a rotation matrix, hence

$\displaystyle\sum_{\chi}a_{\mu\chi}(t)a_{\nu\chi}(t)=\delta_{\mu\nu}$
$\displaystyle\Leftrightarrow\sum\limits_{\chi}\widehat{a}_{\mu\chi}(\Omega)*\widehat{a}_{\nu\chi}(\Omega)=\delta_{\mu\nu}\delta(\Omega)$
(74)

The presented properties of the interaction-frame trajectory hitherto always
hold. However some symmetries of $a_{\mu\nu}(t)$ and
$\widehat{a}_{\mu\nu}(\Omega)$ depend on the basis of the spin operator. In
the hermitian basis $\\{\textrm{I}_{x},\textrm{I}_{y},\textrm{I}_{z}\\}$ the
interaction-frame trajectory is real and therefore the elements of the
frequency-domain interaction-frame trajectory are hermitian functions

$\displaystyle\widehat{a}^{*}_{\mu\nu}(t)$
$\displaystyle=\widehat{a}_{\mu\nu}(t)$ (75)
$\displaystyle\widehat{a}^{*}_{\mu\nu}(-\Omega)$
$\displaystyle=\widehat{a}_{\mu\nu}(\Omega)$ (76)

In contrast in the basis $\\{\textrm{I}_{+},\textrm{I}_{-},\textrm{I}_{z}\\}$
complex conjugation flips the sign of the indices, resulting in

$\displaystyle a^{*}_{\pm\pm}(t)$
$\displaystyle=a_{\mp\mp}(t)\,\quad\quad\textrm{and}\quad
a^{*}_{\pm\mp}(t)\,\,=a_{\mp\pm}(t)$ (77) $\displaystyle
a^{*}_{\pm\pm}(\Omega)$
$\displaystyle=a_{\mp\mp}(-\Omega)\quad\textrm{and}\quad
a^{*}_{\pm\mp}(\Omega)=a_{\mp\pm}(-\Omega)$ (78) $\displaystyle
a^{*}_{z\pm}(t)$
$\displaystyle=a_{z\mp}(t)\quad\textrm{and}\quad\widehat{a}_{z+}(\Omega)=\widehat{a}^{*}_{z-}(-\Omega)$
(79)

####  II.5.2 Repetitive irradiation - nonzero effective field

Many pulse schemes in NMR are repetitive, i.e., consist of repeating segments.
In general a repetitive irradiation does not generate a cyclic interaction-
frame trajectory but will lead to an effective field after each segment. In
this section we consider a pulse scheme of duration $T$, consisting of $N$
repeating segments of duration $\tau$. The interaction-frame trajectory of
such an irradiation can be written as

$\displaystyle a_{\mu\nu}(t)=$ $\displaystyle\begin{dcases}a_{\mu\nu}(t)&0\leq
t\leq\tau\\\ \sum_{\chi}a_{\mu\chi}(t-\tau)a_{\chi\nu}(\tau)&\tau\leq t\leq
2\tau\\\ \sum_{\chi}a_{\mu\chi}(t-2\tau)a_{\chi\nu}(2\tau)&2\tau\leq t\leq
3\tau\\\ \hskip 56.9055pt\vdots&\quad\quad\vdots\\\
\sum_{\chi}a_{\mu\chi}(t-(N-1)\tau)a_{\chi\nu}((N-1)\tau)&(N-1)\tau\leq t\leq
N\tau\end{dcases}$ (80)

Notice that $a_{\mu\nu}(N\tau)=[\mathbf{a}^{N}(\tau)]_{\mu\nu}$ where
$\mathbf{a}(\tau)$ is the matrix with the elements $a_{\mu\nu}(t)$. The
rotation matrix $\mathbf{a}(\tau)$ represents the rotation of the spin due to
the effective field. As shown in detail in the SI, we obtain for a frequency-
domain interaction-frame trajectory

$\displaystyle\widehat{a}_{\mu\nu}(\Omega)$ $\displaystyle=\frac{1}{2\pi
N\tau}\sum_{n=1}^{N}\sum_{\chi}\int_{(n-1)\tau}^{n\tau}a_{\mu\chi}(t-(n-1)\tau)$
$\displaystyle\hskip 85.35826pt\times a_{\chi\nu}((n-1)\tau)e^{-i\Omega t}dt$
(81)
$\displaystyle=\sum_{\chi}\bar{a}_{\mu\chi}(\Omega)\,g^{(N)}_{\chi\nu}(\Omega)$
(82)

with

$\displaystyle
g^{(N)}_{\mu\nu}(\Omega)\coloneqq\frac{1}{N}\sum_{n=0}^{N-1}e^{-in\Omega\tau}a_{\mu\nu}(n\,\tau)$
(83)

and

$\displaystyle\bar{a}_{\mu\nu}(\Omega)$
$\displaystyle=\frac{1}{2\pi\tau}\int\limits_{0}^{\tau}{a}_{\mu\nu}(t)e^{-i\Omega
t}dt$ (84)
$\displaystyle=\mathcal{F}\\{\widetilde{a}_{\mu\nu}(t)\Pi\left(t/\tau-1/2\right)\\}$
(85)
$\displaystyle=\mathcal{F}\\{\widetilde{a}_{\mu\nu}(t)\\}*\mathcal{F}\left\\{\Pi(t/\tau-1/2)\right\\}$
(86)

As a result the $\widehat{a}_{\mu\nu}(\Omega)$ can be calculated efficiently
using DFT, since $a_{\mu\nu}(t)$ can be reduced to the periodized trajectory
of a single segment.

####  II.5.3 Cyclic interaction-frame trajectory - zero effective field

As shown previously, a cyclic interaction-frame trajectory $\mathbf{a}(t)$ can
be rewritten using modulo and a rectangular window function as

$\displaystyle{a}_{\mu\nu}(t)=a_{\mu\nu}(t\,\textrm{mod}\,\tau)\,\Pi\left(\frac{t-T/2}{T}\right)$
(87)

The effective field for a cyclic interaction-frame trajectory after a segment
is always zero, i.e., the interaction frame ends at the same point where it
started. As a consequence $\mathbf{a}(N\tau)=\mathbf{1}$, which leads to

$\displaystyle
g^{(N)}_{\mu\nu}(\Omega)=\frac{1}{N}\sum_{n=0}^{N-1}e^{-in\Omega\tau}$ (88)

Notice, that in the case of a periodic interaction-frame trajectory, i.e., an
infinite repetition of the basic pulse scheme, the relation in Eq. (69)
connects the standard Floquet theory with continuous Floquet theory.

####  II.5.4 Single irradiation axis

In many pulse schemes in NMR, rf irradiation is always along the same axis.
Therefore, the axis of rotation does not change during the pulse scheme if we
do not include the chemical-shift offset into the interaction-frame
calculation. In this case, the phase of the rf-field Hamiltonian can be kept
constant whereas the amplitude $\omega_{1}$ can be positive or negative and
allows rotations in both directions. A special feature of the resulting rf-
field Hamiltonian is that it commutes with itself at different times and,
therefore, an analytical expression for the frequency trajectory
$a_{\mu\nu}(\Omega)$ can be found. Without loss of generality, we consider the
case where $\phi=0$ and, therefore, ${\mathcal{H}}_{\textrm{rf}}$ only
generates rotation around the x-axis. It is convenient to transform in a
tilted frame, where the z-axis is aligned with the rotation axis of the rf-
field. The propagator resulting from ${\mathcal{H}}_{\textrm{rf}}$ together
with the tilted-frame transformation is

$\displaystyle\textrm{U}(t)=e^{-i\int\limits_{0}^{t}\omega_{1}(t)\textrm{d}t\,\textrm{I}_{z}}e^{i\frac{\pi}{2}\textrm{I}_{y}}$
(89)

In most cases we only have to calculate the evolution of the spin operator
along the static magnetic field, which typically is along the z-axis

$\displaystyle\textrm{I}_{z}(t)$
$\displaystyle=\textrm{U}^{\dagger}(t)\textrm{I}_{z}\textrm{U}(t)$
$\displaystyle=\frac{1}{2}\left(e^{i\int\limits_{0}^{t}\omega_{1}(t^{\prime})dt^{\prime}}\,\textrm{I}^{+}+e^{-i\int\limits_{0}^{t}\omega_{1}(t^{\prime})dt^{\prime}}\,\,\textrm{I}^{-}\right)$
$\displaystyle=\frac{1}{2}\left(a_{z+}(t)\,\textrm{I}^{+}+a_{z-}(t)\,\textrm{I}^{-}\right)$
$\displaystyle=\frac{\tau}{2}\int\limits_{-\infty}^{\infty}\left(a_{z+}(\Omega)\,\textrm{I}^{+}+\widehat{a}_{z-}(\Omega)\,\textrm{I}^{-}\right)e^{i\Omega
t}\textrm{d}\Omega$ (90)

As a result, we get a relation between $\omega_{1}(t)$ and
$\widehat{a}_{z\pm}(\Omega)$

$\displaystyle a_{z\pm}(t)=e^{\pm
i\int\limits_{0}^{t}\omega_{1}(t^{\prime})dt^{\prime}}=\tau\int\limits_{-\infty}^{\infty}\widehat{a}_{z\pm}(\Omega)e^{i\Omega
t}\,d\Omega$ (91)
$\displaystyle\widehat{a}_{z\pm}(\Omega)=\frac{1}{2\pi\tau}\int\limits_{-\infty}^{\infty}e^{\pm
i\int\limits_{0}^{t}\omega_{1}(t^{\prime})dt^{\prime}}e^{-i\Omega t}\,dt$ (92)

We can solve Eq. (91) to obtain a solution for $\omega_{1}(t)$

$\displaystyle\omega_{1}(t)$ $\displaystyle=\mp
i\,\frac{d}{dt}\ln\left(\mathcal{F}^{-1}\\{\widehat{a}_{z\pm}(\Omega)\\}\right)$
$\displaystyle=\pm\frac{\mathcal{F}^{-1}\\{\Omega\,\widehat{a}_{z\pm}(\Omega)\\}}{\mathcal{F}^{-1}\\{\widehat{a}_{z\pm}(\Omega)\\}}$
(93)

In addition, from the convolution theorem we obtain

$\displaystyle\widehat{a}_{\mu\nu}^{(1)}(\Omega)*\widehat{a}_{\mu\nu}^{(2)}(\Omega)*\ldots*\widehat{a}_{\mu\nu}^{(n)}(\Omega)$
$\displaystyle=\mathcal{F}\biggr{\\{}\prod\limits_{j=1}^{n}\,a_{\mu\nu}^{(j)}(t)\biggl{\\}}$
$\displaystyle=\mathcal{F}\biggr{\\{}\exp\biggr{(}$
$\displaystyle-i\,\int_{0}^{t}\sum\limits_{j=1}^{n}\omega_{1}^{(j)}(t)\textrm{d}t\biggl{)}\biggl{\\}}$
(94)

This means that, if we decompose the desired shape
$\widehat{a}_{\mu\nu}(\Omega)$ as a convolution, we can just add the
corresponding $\omega_{1}(t)$ to obtain the desired rf-field profile.

Figure 1: Schematic diagram of continuous Floquet theory for magnetic
resonance. The information of the interaction-frame transformation
$\textrm{U}(t)$ is encoded in the frequency-domain interaction-frame
trajectory $\widehat{\mathbf{a}}(\Omega)$ which is calculated from the
interaction-frame trajectory $\mathbf{a}(t)$. Subsequently
$\widehat{\mathbf{a}}(\Omega)$ is combined with the effective Hamiltonian.

###  II.6 Calculation of the effective-field Hamiltonian from the frequency-
domain interaction-frame trajectory

In this section we reverse the previous procedure and calculate the rf-field
Hamiltonian from $\widehat{a}_{\mu\nu}(\Omega)$. As a first step we simply
apply the inverse Fourier transform

$\displaystyle
a_{\mu\nu}(t)=\tau\int\limits_{-\infty}^{\infty}\widehat{a}_{\mu\nu}(\Omega)e^{i\Omega
t}\,d\Omega$ (95)

In the next step we want to construct $\textrm{U}(t)$ from ${\bf{a}}(t)$. As
mentioned previously the U(t) given in Eq. (55) are members of $SU(2)/Z_{2}$,
which is isomorphic to $SO(3)$. Hence we can map each rotation matrix
${\bf{a}}_{j}$ uniquely to $\textrm{U}_{j}$. For the construction, we use
quaternions, Cartan (1966); Blümich and Spiess (1985) since they can represent
complex as well as real rotations. More precisely, we reformulate the rotation
matrices as unit quaternions and subsequently represent them with spin
matrices in order to read out the pulse parameter.

A unit quaternion is given by $q_{j}=u+j\,v+k\,w+l\,z$ with $j\,k\,l=-1$ and
$u^{2}+v^{2}+w^{2}+z^{2}=1$. With the Euler-Rodriguez formula, Rodrigues
(1840); Cartan (1966) we express its components with the elements of the
corresponding rotation matrix

$\displaystyle u$ $\displaystyle=\frac{1}{2}\sqrt{1+tr({\mathbf{a}})}$ (96)
$\displaystyle v$ $\displaystyle=\frac{1}{4\,a}(a_{zy}-a_{yz})$ (97)
$\displaystyle w$ $\displaystyle=\frac{1}{4\,a}(a_{xz}-a_{zx})$ (98)
$\displaystyle z$ $\displaystyle=\frac{1}{4\,a}(a_{yx}-a_{xy})$ (99)

Notice, that there are different options to calculate a unit quaternion from
the corresponding rotation matrix. It is important for the numerical
evaluation to choose the option, where the denominator is not close to zero,
in order to increase the precision. We can always represent a unit quaternion
as

$\displaystyle q_{j}=\cos\left(\frac{|\theta|\delta
t}{2}\right)+\frac{1}{|\theta|}(j\,\theta_{x}+k\,\theta_{y}+l\,\theta_{z})\sin\left(\frac{|\theta|\delta
t}{2}\right)$ (100)

where $\vec{\theta}(t)$ is defined in Eq. (53). The norm of $\vec{\theta}(t)$
is the angular velocity and the direction specifies the rotation axis. Next,
we rewrite the unit quaternions to assemble $\textrm{U}_{j}$ in its
exponential form. Therefore, we represent the quaternion units with the spin
matrices

$\displaystyle j=2i\,\textrm{I}_{x},\hskip 8.5359ptk=2i\,\textrm{I}_{y},\hskip
8.5359ptl=2i\,\textrm{I}_{z}$ (101)

which leads to the desired representation

$\displaystyle q_{j}$
$\displaystyle=\mathbb{I}\,\cos\left(\frac{|\theta|\delta
t}{2}\right)-i\frac{\vec{\theta}}{|\theta|}2\,\vec{\textrm{I}}\,\sin\left(\frac{|\theta|\delta
t}{2}\right)$ $\displaystyle=\exp(-i\,\vec{\theta}\,\,\vec{\textrm{I}}\,\delta
t)=\exp(-i\,{\mathcal{H}}_{\textrm{rf}}\,\delta t)$ (102)

Together with the definition of $\vec{\theta}$ in Eq. (53) we obtain

$\displaystyle u$
$\displaystyle=\frac{1}{2}\sqrt{1+\textrm{tr}(a)}=\cos\left(\frac{|\theta|\delta
t}{2}\right)$ (103) $\displaystyle v$
$\displaystyle=\frac{1}{4\,u}(a_{zy}-a_{yz})=\frac{\omega_{1}}{|\theta|}\cos(\phi)\,\sin\left(\frac{|\theta|\delta
t}{2}\right)$ (104) $\displaystyle w$
$\displaystyle=\frac{1}{4\,u}(a_{xz}-a_{zx})=\frac{\omega_{1}}{|\theta|}\sin(\phi)\,\sin\left(\frac{|\theta|\delta
t}{2}\right)$ (105) $\displaystyle z$
$\displaystyle=\frac{1}{4\,u}(a_{yx}-a_{xy})=\frac{\omega_{z}}{|\theta|}\,\sin\left(\frac{|\theta|\delta
t}{2}\right)$ (106)

The solution of these equations for the pulse parameter
$\omega_{1},\,\phi,\,\omega_{z}$ are

$\displaystyle|\theta|\delta t$ $\displaystyle=2\,\arccos\left(u\right)$ (107)
$\displaystyle\omega_{1}$
$\displaystyle=\frac{|\theta|\,w}{\sin(\phi)\,\sin\left(\frac{|\theta|\delta
t}{2}\right)}$ (108) $\displaystyle\phi$
$\displaystyle=\arctan\left(\frac{w}{v}\right)$ (109)
$\displaystyle\omega_{z}$
$\displaystyle=\frac{|\theta|\,z}{\sin\left(\frac{|\theta|\delta
t}{2}\right)}$ (110)

which can be also expressed using $a_{\mu\nu}$ as

$\displaystyle|\theta|\delta
t=\arccos\left(\frac{1}{2}\sqrt{1+\textrm{tr}({\bf{a}})}\right)$ (111)
$\displaystyle\phi=\arctan\left(\frac{a_{xz}-a_{zx}}{a_{zy}-a_{yz}}\right)$
(112)
$\displaystyle\omega_{1}=-\frac{\text{arcsec}\left(2\,(\textrm{tr}({\bf{a}})+1)^{-\frac{1}{2}}\right)(a_{xz}-a_{zx}+a_{yz}-a_{zy})^{2}\csc\left(2\,\text{arcsec}\left(2\,(\textrm{tr}({\bf{a}})+1)^{-\frac{1}{2}}\right)\right)}{\delta
t\left[\left(\textrm{tr}({\bf{a}})+1\right)\left((a_{xz}-a_{zx})^{2}+(a_{yz}-a_{zy})^{2}\right)\right]^{\frac{1}{2}}}$
(113)
$\displaystyle\omega_{z}=\frac{2\,(a_{yx}-a_{xy})\,\text{arcsec}\left(2\,(\textrm{tr}({\bf{a}})+1)^{-\frac{1}{2}}\right)}{\delta
t\left[\left(\textrm{tr}(\bf{a})+1\right)\left(3-\textrm{tr}({\bf{a}})\right)\right]^{\frac{1}{2}}}$
(114)

Eq. (112 \- 114) express the phase, rf-field amplitude and offset as function
of the interaction frame trajectory $a_{\mu\nu}(t)$. In order to obtain the
complete sequence, these equations have to be evaluated at every time step of
$a_{\mu\nu}(t_{j})$.

Figure 1 depicts the procedure to obtain the frequency-domain interaction-
frame trajectory from the pulse parameter and vice versa. The mapping between
the pulse scheme parameter frequency-domain interaction-frame trajectory is
bijective, hence we always can find an unique result in both directions.
However, only if the conditions given by Eqs. (77-79) are fulfilled by the
frequency-domain interaction-frame trajectory, a physical solution for the
pulse parameter can be found.

##  III The design of band-selective MIRROR experiments

In this section we apply continuous Floquet theory to design the rf
irradiation for a tunable PDSD-based sequence called AM-MIRROR. Scholz  _et
al._ (2008); Wittmann  _et al._ (2014, 2016) The sequence can be designed for
broadband or band-selective polarization transfer even at fast MAS frequencies
by creating side-bands in the zero-quantum spectrum. Using the continuous
frequency-space Floquet formalism, we show how to calculate the irradiation
scheme to achieve tailored zero-quantum side-bands without any need of
numerical optimization, which is in general not possible with the standard
Floquet approach.

###  III.1 Theoretical preliminaries

To analyze the MIRROR experiment it is sufficient to consider a
$\textrm{S}_{2}\textrm{I}$ spin system with irradiation on the I spin. In the
standard rotating frame such a Hamiltonian has the form

$\displaystyle{\mathcal{H}}(t)={\mathcal{H}}_{S}(t)+{\mathcal{H}}_{\textrm{rf}}(t)$
(115)

with

$\displaystyle{\mathcal{H}}_{S}(t)=$
$\displaystyle\sum\limits_{n=-2}^{2}\omega_{S_{1}S_{2}}^{(n)}e^{in\omega_{\textrm{r}}t}\,(3\,\textrm{S}_{1z}\textrm{S}_{2z}-\vec{\textrm{S}}_{1}\vec{\textrm{S}}_{2})$
$\displaystyle+\sum\limits_{p=1}^{2}\sum\limits_{n=-2}^{2}\omega_{I_{1}S_{p}}^{(n)}e^{in\omega_{\textrm{r}}t}\,2\,\textrm{I}_{z}\textrm{S}_{pz}$
$\displaystyle+\sum\limits_{p=1}^{2}\omega_{S_{p}}^{(0)}\textrm{S}_{pz}$ (116)

and

$\displaystyle{\mathcal{H}}_{\textrm{rf}}(t)=\omega_{1}(t)\,\textrm{I}_{x}$
(117)

Here, consider only isotropic the chemical shifts of the S spins, but not of
the I spins. Furthermore we assume an amplitude-modulated rf-irradiation in
the x-direction ($\phi=0,180^{\circ}$). For simplicity, we tilt the frame of
reference to align the z-axis with rf-irradiation axis as in section II.5.4.
In addition, we transform into an interaction frame with the isotropic
chemical shift of the S spins. As usual we also transform into an interaction
frame with the rf irradiation on the I spins leading to a total interaction-
frame transformation defined by:

$\displaystyle\textrm{U}(t)=$
$\displaystyle\underbrace{{\color[rgb]{0,0,0}\exp\biggl{(}-i\int\limits_{0}^{t}\omega_{1}(t^{\prime})dt^{\prime}\,{\color[rgb]{0,0,0}\textrm{I}_{z}}\biggl{)}}}_{\textrm{rf-
irradiation}}{\color[rgb]{0,0,0}\underbrace{\vphantom{\int\limits_{0}^{t}}{\color[rgb]{0,0,0}\exp\biggl{(}-i\frac{\pi}{2}\,\textrm{I}_{y}\biggl{)}}}_{\textrm{tilted-
frame}}}$ (118)
$\displaystyle\underbrace{\vphantom{\int\limits_{0}^{t}}{\color[rgb]{0,0,0}\exp{\biggl{(}i\sum\limits_{p=1}^{2}\omega_{Sp}^{(0)}t\,\textrm{S}_{pz}\biggl{)}}}}_{\textrm{chemical
shift}}$

Since the S-spin chemical shift and the MAS spinning lead to a continuous
rotation around a single axis, the interaction-frame Hamiltonian
$\widetilde{{\mathcal{H}}}(t)\coloneqq\textrm{U}(t){\mathcal{H}}(t)\textrm{U}^{\dagger}(t)$
can be written as:

$\displaystyle\widetilde{{\mathcal{H}}}(t)=T\sum\limits_{n=-2}^{2}\sum\limits_{\ell=-1}^{1}\,\int\limits_{-\infty}^{\infty}{\mathcal{H}}^{(n,\ell)}(\Omega)\,e^{in\omega_{\textrm{r}}t}e^{i\ell\Delta\omega_{\textrm{iso}}t}e^{i\Omega
t}\,\textrm{d}\Omega$ (119)

from where we can read out the MIRROR resonance condition

$\displaystyle\Omega_{0}+n_{0}\omega_{r}+\ell_{0}\Delta\omega_{\textrm{iso}}=0$
(120)
$\displaystyle(\Omega_{0},n_{0},\ell_{0})\in\mathbb{R}\times\\{\textrm{-}2,\textrm{-}1,...,2\\}\times\\{\textrm{-}1,0,1\\}$
(121)

The dominant second-order effective Hamiltonian is

$\displaystyle\widetilde{{\mathcal{H}}}_{SS\otimes
IS}=\frac{1}{4}\sum\limits_{n,\nu=-2}^{2}\,\int\limits_{-\infty}^{\infty}$
$\displaystyle\biggl{(}\frac{\omega_{S_{1}S_{2}}^{(\nu)}(\omega_{I_{1}S_{2}}^{(n-\nu)}-\omega_{I_{1}S_{2}}^{(n-\nu)})}{\Delta\omega_{\textrm{iso}}-\nu\omega_{\textrm{r}}}+\frac{\omega_{S_{1}S_{2}}^{(n-\nu)}(\omega_{I_{1}S_{2}}^{(\nu)}-\omega_{I_{1}S_{2}}^{(\nu)})}{\Delta\omega_{\textrm{iso}}+\nu\omega_{\textrm{r}}}\biggr{)}\,(\widehat{a}_{z\pm}(\Omega)\,\textrm{I}^{\pm}+\widehat{a}_{z\mp}(\Omega)\,\textrm{I}^{\mp})\,\textrm{S}_{1}^{\pm}\textrm{S}_{2}^{\mp}\,d\Omega$
(122)

where

$\displaystyle\widehat{a}_{z\pm}(\Omega)$
$\displaystyle\coloneqq\widehat{a}_{zx}(\Omega)\pm
i\,\widehat{a}_{zy}(\Omega)$ (123)

In the following we are going to use Eq. (93) which enables the calculation of
the rf-field amplitude.

###  III.2 Targeting single $\Delta\omega_{\textrm{iso}}$

The first example is the recoupling of a single chemical-shift difference
$\Delta\omega_{\mathrm{iso}}$ given by

$\displaystyle\widehat{a}_{z+}(\Omega)\stackrel{{\scriptstyle!}}{{=}}\delta(\Omega-\Delta\omega_{\textrm{iso}})\hskip
20.00003pt$ (124)

Inserting it in Eq. (91) leads to

$\displaystyle
e^{i\int\limits_{0}^{t}\omega_{1}(t^{\prime})dt^{\prime}}=\int\limits_{-\infty}^{\infty}\delta(\Omega-\Delta\omega_{\textrm{iso}})e^{i\Omega
t}\,d\Omega=e^{i\Delta\omega_{\textrm{iso}}t}$ (125)

hence

$\displaystyle\omega_{1}(t)=\Delta\omega_{\textrm{iso}}$ (126)

Alternatively we can use Eq. (93)

$\displaystyle\omega_{1}(t)$
$\displaystyle=\frac{\mathcal{F}^{-1}\\{\Omega\,\delta(\Omega-\Delta\omega_{\textrm{iso}})\\}}{\mathcal{F}^{-1}\\{\delta(\Omega-\Delta\omega_{\textrm{iso}})\\}}=\frac{\Delta\omega_{\textrm{iso}}e^{-i\Delta\omega_{\textrm{iso}}t}}{e^{-i\Delta\omega_{\textrm{iso}}t}}=\Delta\omega_{\textrm{iso}}$
(127)

This solution is not physical, since it would require irradiation of a
infinite duration. As we will show below, the optimal rf-field amplitude for a
given duration $T$ is in our case

$\displaystyle\omega_{1}(t)=\begin{cases}\Delta\omega_{\textrm{iso}}&\,0\leq
t\leq T\\\ 0&\hskip 15.00002ptelse\end{cases}$ (128)

Using Eq. (92), we obtain

$\displaystyle\widehat{a}_{z+}(\Omega)$ $\displaystyle=\frac{i}{2\pi
T}\,\frac{1-e^{i\text{$T$}(\Omega-\text{$\Delta\omega_{\textrm{iso}}$})}}{\Delta\omega_{\textrm{iso}}-\Omega}$
(129)

Note that this is a sinc function where the time origin was shifted to the
beginning of the time period and not at the center. The singularity at
$\Omega=\Delta\omega_{\textrm{iso}}$ is removable, which can be seen from its
expansion

$\displaystyle\widehat{a}_{z+}(\Omega)$
$\displaystyle=\frac{i}{2\pi}\,\sum\limits_{m=1}^{\infty}\frac{i^{m}}{m!}[T(\Delta\omega_{\textrm{iso}}-\Omega)]^{m-1}$
(130)

We can easily calculate the global maximum of $|\widehat{a}_{z+}(\Omega)|$
with Eq. (130), which is located at $\Omega=\Delta\omega_{\textrm{iso}}$

$\displaystyle|a_{z+}(\Delta\omega_{\textrm{iso}})|$
$\displaystyle=\frac{1}{2\pi}$ (131)

Notice, that the global maximum is independent of the duration of the pulse
scheme, because of the chosen normalization.

###  III.3 Targeting multiple isolated $\Delta\omega_{\textrm{iso}}$

The simplest example of targeting multiple frequencies is to target 2
frequencies

$\displaystyle\widehat{a}_{z+}(\Omega)\stackrel{{\scriptstyle!}}{{=}}\frac{1}{2}\left(\delta(\Omega-\Delta\omega_{\textrm{iso}}^{(1)})+\delta(\Omega-\Delta\omega_{\textrm{iso}}^{(2)})\right)$
(132)

Inserting it into Eq. (93) leads to

$\displaystyle\omega_{1}(t)=\Delta\omega_{\textrm{iso}}^{(1)}+\Delta\omega_{\textrm{iso}}^{(2)}$
(133)

This irradiation is not possible for 2 reasons. First, we would again require
irradiation of infinite duration as before. The second problem is that two
waves with the same phase just add up in amplitude, i.e.
$\Delta\omega_{\textrm{iso}}^{(1)}+\Delta\omega_{\textrm{iso}}^{(2)}=\Delta\omega_{\textrm{iso}}^{(3)}$.
As a result, we would just target one frequency, namely
$\Delta\omega_{\textrm{iso}}^{(3)}$. Notice, that this is the only existing
solution for the given $\widehat{a}_{z+}(\Omega)$. Therefore we have no choice
than to split the irradiation into parts that target the two frequencies
separately and successively in time

$\displaystyle\omega_{1}(t)=\begin{cases}\Delta\omega_{\textrm{iso}}^{(1)}&\,0\leq
t\leq\tau_{1}\\\ \Delta\omega_{\textrm{iso}}^{(2)}&\tau_{1}<t\leq
T\end{cases}$ (134)

We again can use Eq . (92) and obtain

$\displaystyle\widehat{a}_{z+}(\Omega)$
$\displaystyle=\widehat{a}^{(1)}_{z+}(\Omega)+\widehat{a}^{(2)}_{z+}(\Omega)$
$\displaystyle=\frac{i}{2\pi
T}\frac{\left(1-e^{-i\text{$\tau_{1}$}(\Omega-\text{$\Delta\omega_{\textrm{iso}}^{(1)}$})}\right)}{(\text{$\Delta\omega_{\textrm{iso}}^{(1)}$}-\Omega)}$
$\displaystyle\quad+\frac{i}{2\pi
T}\frac{\left(e^{-i\text{$\tau_{1}$}(\Omega-\text{$\Delta\omega_{\textrm{iso}}^{(2)}$})}-e^{-i\text{$T$}(\Omega-\text{$\Delta\omega_{\textrm{iso}}^{(2)}$})}\right)e^{i(\text{$\tau_{1}$}\text{$\Delta\omega_{\textrm{iso}}^{(1)}$}-\text{$T$}\text{$\Delta\omega_{\textrm{iso}}^{(2)}$})}}{(\text{$\Delta\omega_{\textrm{iso}}^{(2)}$}-\Omega)}$
(135)

Notice, that only the phase factor of the second term $a^{(2)}_{+}(\Omega)$
depends on the first part of the irradiation, but not its absolute value.
Similar to the case of a single chemical shift we find

$\displaystyle|\widehat{a}^{(1)}_{z+}(\Delta\omega^{(1)}_{\textrm{iso}})|$
$\displaystyle=\frac{\tau_{1}}{2\pi T}$ (136)
$\displaystyle|\widehat{a}^{(2)}_{z+}(\Delta\omega^{(2)}_{\textrm{iso}})|$
$\displaystyle=\frac{T-\tau_{1}}{2\pi T}$ (137)

which approximates well the height of the two highest maxima of
$|\widehat{a}_{z+}(\Omega)|$, becoming more accurate the bigger
$|\Delta\omega^{(2)}_{\textrm{iso}}-\Delta\omega^{(1)}_{\textrm{iso}}|$,
because the overlap between $a^{(1)}(\Omega)$ and $a^{(2)}(\Omega)$ decreases.
We can generalize this findings to $N$ chemical shifts as

$\displaystyle\omega_{1}(t)=\begin{dcases}\begin{aligned}
&\Delta\omega^{(1)}_{\textrm{iso}}\hskip 23.99997pt0\leq t\leq\tau_{1}\\\
&\Delta\omega^{(2)}_{\textrm{iso}}\hskip
20.00003pt\tau_{1}<t\leq\tau_{2}\\\\[-4.30554pt] &\hskip
10.00002pt\mathmakebox[\widthof{{}\Omega{}}][c]{\vdots}\hskip
50.00008pt\vdots\\\\[-2.15277pt] &\Delta\omega^{(n)}_{\textrm{iso}}\hskip
8.00003pt\tau_{n-1}<t\leq\tau_{n}\end{aligned}\end{dcases}$ (138)

leading to

$\displaystyle\widehat{a}_{z+}(\Omega)$ $\displaystyle=\frac{i}{2\pi
T}\sum\limits_{n=0}^{N}\frac{\left(e^{-i\text{$\tau_{n-1}$}(\Omega-\text{$\omega_{n}$})}-e^{-i\tau_{n}(\Omega-\text{$\omega_{n}$})}\right)}{(\omega_{n}-\Omega)}$
$\displaystyle\quad
e^{i(\tau_{1}\omega_{1}+\tau_{2}\omega_{2}+\cdots+\tau_{n-2}\omega_{n-2}+\tau_{n-1}(\omega_{n-1}-\omega_{n}))}$
(139)

where we set $\tau_{N}=T$ and $\tau_{0}=0=\omega_{0}$. Like before, the maxima
of the absolute value of each term is given by

$\displaystyle|\widehat{a}^{(n)}_{z+}(\Delta\omega^{(n)}_{\textrm{iso}})|$
$\displaystyle=\frac{\tau_{n+1}-\tau_{n}}{2\pi T}$ (140)

It might be possible to construct an arbitrary function from general stepwise
function by taking the limits carefully. The same limit can be applied
directly on Eq. (139), if the integral in Eq. (91) is dominant integrable
Bartle (2014) and hence we can change the order of the limit and the intergral
operation.

###  III.4 Targeting a range of $\Delta\omega_{\textrm{iso}}$

To obtain $\omega_{1}(t)$ we set $\widehat{a}_{z+}(\Omega)$ to be a
rectangular function with width $b$, normalized according to Eq. (67), as

$\displaystyle\widehat{a}_{z+}(\Omega)\stackrel{{\scriptstyle!}}{{=}}\sqrt{\frac{2\pi}{b}}\,\Pi(\Omega/b)$
(141)

From Eq. (93) follows

$\displaystyle\omega_{1}(t)=\frac{\mathcal{F}^{-1}\\{\Omega\,\,\Pi(\Omega/b)\\}}{\mathcal{F}^{-1}\\{\,\Pi(\Omega/b)\\}}=i\left(\frac{b}{2}\cot\left(\frac{b\,t}{2}\right)-\frac{1}{t}\right)$
(142)

Notice that $\omega_{1}(t)$ is independent of a constant factor of
$\widehat{a}_{z+}(\Omega)$. As before we have the problem of infinite duration
in order to obtain the desired shape exactly.
In contrast to the previous theoretical rf-field amplitudes, we get a time-
dependent function, which has poles and is purely imaginary. Motivated by Eq.
(138) and Eq. (139), we use a linear function

$\displaystyle\omega_{1}(t)=\frac{(\Delta\omega_{\textrm{iso}}^{(2)}-\Delta\omega_{\textrm{iso}}^{(1)})}{\vphantom{T^{5}}T}\,t+\Delta\omega_{\textrm{iso}}^{(1)}$
(143)

Inserting it in Eq. (92), yields

$\displaystyle\widehat{a}_{z+}(\Omega)$
$\displaystyle=\sqrt{\frac{\pi}{2T}}\,e^{i\frac{3\pi}{4}}e^{i\frac{T(\Delta\omega^{(1)}_{\textrm{iso}}-\Omega)^{2}}{2(\Delta\omega^{(1)}_{\textrm{iso}}-\Delta\omega^{(2)}_{\textrm{iso}})}}$
$\displaystyle\times$
$\displaystyle\frac{\text{erf}\left(\frac{\left(\frac{1}{2}+\frac{i}{2}\right)\sqrt{T}(\Delta\omega^{(2)}_{\textrm{iso}}-\Omega)}{\sqrt{\Delta\omega^{(1)}_{\textrm{iso}}-\Delta\omega^{(2)}_{\textrm{iso}}}}\right)-\text{erf}\left(\frac{\left(\frac{1}{2}+\frac{i}{2}\right)\sqrt{T}(\Delta\omega^{(1)}_{\textrm{iso}}-\Omega)}{\sqrt{\Delta\omega^{(1)}_{\textrm{iso}}-\Delta\omega^{(2)}_{\textrm{iso}}}}\right)}{\sqrt{\Delta\omega^{(1)}_{\textrm{iso}}-\Delta\omega^{(2)}_{\textrm{iso}}}}$
(144)

Taking a closer look on Eq. (143) and Eq. (144), we notice that Eq. (143)
consists of 2 parts. The first term causes a rectangular like shape, the
second term shifts this shape to the desired position. An intuitive picture
offers Eq. (91) when we consider an interaction-frame rotation with an
uniformly changing frequency $a_{z+}(t)=e^{i(m\,t+c)t}$. The frequencies
$\omega(t)=m\,t+c$ are visible in its Fourier transform. Sweeping uniformly
through a range of nutation frequencies slowly once, will drive more
resonances than to sweep the same range several times, even if the overall
time is the same. This is because a repeated sweep, also includes other, much
slower, frequency components, due to the repetition.

Figure 2 shows
$\frac{1}{4\pi}\bigl{(}|\widehat{a}_{z+}(\Omega)|^{2}+|\widehat{a}_{z-}(\Omega)|^{2}\bigr{)}$
for different pulse schemes and different number of repetitions of the pulse
scheme. For Figs. (2A) and (2B) a step-wise pulse schemes (Eq. (134) and Eq.
(135)) was used with $\tau_{1}=T/2$ for (A) and $\tau_{1}=3/4T$ for (B). For
Fig. (2C) a ramped irradiation described with Eqs. (143) and (144) was used.
For the red curves the pulse scheme was applied only once during $T=20\,$ms
leading to 4 individual peaks for (A) and (B) and a box-shaped curve for (C).
For the blue curves the irradiation scheme has been applied 10 times during
the same duration of $T=20\,$ms, introducing the modulation frequency
$\tau_{m}=2\,$ms. As a consequence the blue curves is discretized due to the
modulation frequency $\nu_{m}=1/\tau_{m}=0.5\,$kHz. The green curves is
obtained by applying the irradiation scheme once for $T=2\,$ms and scaling the
curves by a factor of 10, which leads to the envelope of the blue curves.
Notice that all the curves shown were obtained by analytical expressions, i.e
using Eq. (135) for (A) and (B) and Eq. (134) for (C). To account for the
effect of the repetition Eq. (83) was used, and therefore all plots are
obtained by FFT of a single element.

Figure 2:
$\frac{1}{4\pi}\bigl{(}|\widehat{a}_{z+}(\Omega)|^{2}+|\widehat{a}_{z-}(\Omega)|^{2}\bigr{)}$
for different irradiation schemes. (A) Step-wise irradiation (Eqs. (134) and
(135)) with $\Delta\omega^{(1)}_{\textrm{iso}}=5$ and
$\Delta\omega^{(2)}_{\textrm{iso}}=10$ with $T=2\tau_{1}=20\,$ms. (B) Same as
(A) but with $\tau_{1}=3/4\,T$. (C) Ramped irradiation described by Eq. (143)
and Eq. (144), ranging from $7.5\,$kHz $\leq\omega_{1}\leq 12.5\,$kHz. For the
red curves, the duration has been set to $T=20$, repeating the shape once. For
the blue curves the duration has been set to $T/10$ and repeated $10$ times.
For the green curves the duration has been $T/10$ and repeated once, but
subsequently scaled up $10$ times. The integral of the red and blue curve is
equal, since the duration of the irradiation is the same. As a consequence of
the repetition, the blue curve is discretized by the modulation frequency
$\nu_{m}=0.5\,\textrm{kHz}$.

###  III.5 Comparison to numerical simulation

An alternate description of such PDSD-type experiments which does not rely on
effective Hamiltonians, was provided by Veshtort Veshtort and Griffin (2011),
where the rate constants $k_{\mathrm{D}}$ under MAS can be calculated as

$\displaystyle
k_{\mathrm{D}}=\sum^{2}_{n=-2}|\omega^{(n)}_{S_{1}S_{2}}|^{2}\,Re\\{J_{0}(\Delta\omega_{\text{iso}}-n\omega_{r})\\}$
(145)

Here $\omega^{(n)}_{S_{1}S_{2}}$ is defined as in Eq. (III.1) and
$J_{0}(\omega)$ is the ZQ line Veshtort and Griffin (2011). Notice, that
$k_{\mathrm{D}}$ is proportional to $J_{0}(n\omega_{\mathrm{r}})$, which can
be manipulated by the irradiation and is dependent on the MAS frequency,
similarly to $\widehat{a}_{z\mp}(\Omega)$. However the zero-quantum line and
$\widehat{a}_{z\mp}(\Omega)$ are quantities that stem from fundamentally
different theoretical descriptions. Nevertheless both quantities ultimately
describe the evolution of the ZQ operator and dictate the ZQ transfer
selectivity and efficiency. Figure 3 shows a comparison of the ZQ line and
$1/(4\pi)(|\widehat{a}_{z+}(\Omega)|^{2}+|\widehat{a}_{z-}(\Omega)|^{2})$ for
different irradiation schemes. The plots of the ZQ line are the results of
numerical time-slicing simulation of a CH2 spin system using the C++ library
GAMMA Smith  _et al._ (1994). The mixing time was set to $20\,$ms and the MAS
frequency to $50\,$kHz. The $\widehat{a}_{z\pm}(\Omega)$ were calculated from
Eq. (135) and (144) using the same mixing time as in the numerical simulation.
As mentioned before the two quantities stem from different descriptions, but
the agreement, especially of the width and position, is very good, underlying
their physical similarities.

Figure 3: Comparison between numerical simulation of the zero quantum line and
the interaction frame trajectory for different irradiation schemes. (A) Step
irradiation as given in Eq. (134) with
$\Delta\omega_{\text{iso}}^{(1)}=5\,$kHz, $\Delta\omega_{\text{iso}}^{(2)}$
and $\tau_{1}=T/2=10$ ms. (B) Ramp irradiation with $7.5\,$kHz
$\leq\omega_{1}\leq 12.5\,$kHz. (C) Ramp irradiation with $7.5\,$kHz
$\leq\omega_{1}\leq 12.5\,$kHz.

##  IV Conclusion

We presented a generalization of operator-based Floquet theory to non-periodic
Hamiltonians. Instead of a discrete frequency space, which is sufficient for
the description of periodic Hamiltonians as in the standard Floquet approach,
we utilized a continuous frequency space allowing non-periodic Hamiltonians
with finite duration. In contrast to periodic modulations, modulation with
finite duration led to resonance conditions that are not infinitely sharp,
enabling the description of non-resonant behaviors within the same framework,
which is very cumbersome in standard Floquet theory. Since the framework is
not restricted to periodic modulations, it can be use in combination with any
interaction-frame transformation. This might be useful to analyze the effect
of specific terms in the Hamiltonian by choosing an appropriate interaction
frame, or simplifying the description by incorporating all single-spin
modulations in the interaction frame. Despite the finite duration of the
modulations, we can still utilize FFT for computationally efficient
calculations and express the effect of the repetition of a pulse scheme by a
simple quantity. We saw that the pulse parameter and the interaction-frame
trajectory, which encodes the modulation of the spin system, can be mapped to
each other bijectively, enabling reverse engineering pulse schemes from the
effective Hamiltonian to a certain extent. Finally we applied the formalism to
a PDSD-based sequence called AM-MIRROR where we could tailor desired zero-
quantum side-bands using the mentioned bijective mapping, without any need of
numerical optimization. This framework was developed to describe and design
solid state NMR experiments under MAS, but is not restricted to it. It can be
used to describe experiments of other magnetic resonance technologies such as
DNP and EPR.

## ACKNOWLEDGMENTS

This research has been supported by the ETH Zürich and the Schweizerischer
Nationalfonds zur Förderung der Wissenschaftlichen Forschung (grant no.
200020_188988).

## AUTHORS DECLARATION

### Conflict of interest

The authors have no conflicts of interest to disclose.

## DATA AVAILABILITY

The simulated data will be uploaded to a public repository after acceptance of
the paper.

## References

  * Andrew, Bradbury, and Eades (1958) E. R. Andrew, A. Bradbury, and R. G. Eades, “Nuclear Magnetic Resonance Spectra From a Crystal Rotated at High Speed,” Nature 182, 1659–1659 (1958).
  * Andrew, Bradbury, and Eades (1959) E. R. Andrew, A. Bradbury, and R. G. Eades, “Removal of Dipolar Broadening of Nuclear Magnetic Resonance Spectra of Solids by Specimen Rotation,” Nature 183, 1802–1803 (1959).
  * Lowe (1959) I. J. Lowe, “Free Induction Decays of Rotating Solids,” Phys. Rev. Lett. 2, 285–287 (1959).
  * Samoson, Lippmaa, and Pines (1988) A. Samoson, E. Lippmaa, and A. Pines, “High resolution solid-state N.M.R.” Mol. Phys. 65, 1013–1018 (1988).
  * Mueller  _et al._ (1990) K. T. Mueller, B. Sun, G. CHINGAS, J. W. Zwanziger, T. Terao, and A. Pines, “Dynamic-Angle Spinning of Quadrupolar Nuclei,” J. Magn. Reson. 86, 470–487 (1990).
  * Haeberlen (1976) U. Haeberlen,  _High resolution NMR in solids : selective averaging_ (New York : Academic Press, 1976).
  * Mehring (1983) M. Mehring,  _Principles of high resolution NMR in solids_ (Springer-Verlag, 1983).
  * Haeberlen and Waugh (1968) U. Haeberlen and J. S. Waugh, “Coherent Averaging Effects in Magnetic Resonance,” Phys. Rev. 175, 453–467 (1968).
  * Maricq (1982) M. M. Maricq, “Application of Average Hamiltonian Theory to the Nmr of Solids,” Phys. Rev. B 25, 6622–6632 (1982).
  * Ernst, Bodenhausen, and Wokaun (1990) R. R. Ernst, G. Bodenhausen, and A. Wokaun,  _Principles of Nuclear Magnetic Resonance in One and Two Dimensions_ (Oxford University Press, Oxford, 1990).
  * Floquet (1883) G. Floquet, “Sur la théorie des équations différentielles linéaires,” Annales Scientifiques de l’E.N.S. 12, 47–88 (1883).
  * Shirley (1965) J. H. Shirley, “Solution of Schrodinger Equation with a Hamiltonian Periodic in Time,” Phys. Rev. 138, B979–B987 (1965).
  * Scholz, van Beek, and Ernst (2010) I. Scholz, J. D. van Beek, and M. Ernst, “Operator-based Floquet theory in solid-state NMR,” Solid State Nucl. Magn. Reson. 37, 39–59 (2010).
  * Leskes, Madhu, and Vega (2010) M. Leskes, P. K. Madhu, and S. Vega, “Floquet theory in solid-state nuclear magnetic resonance.” Prog. NMR Spectr. 57, 345–380 (2010).
  * Ivanov  _et al._ (2021) K. L. Ivanov, K. R. Mote, M. Ernst, A. Equbal, and P. K. Madhu, “Floquet theory in magnetic resonance: Formalism and applications,” Prog. NMR Spectr. 126-127, 17–58 (2021).
  * Madhu and Kurur (2006) P. K. Madhu and N. D. Kurur, “Fer expansion for effective propagators and Hamiltonians in NMR,” Chem. Phys. Lett. 418, 235–238 (2006).
  * Takegoshi  _et al._ (2015) K. Takegoshi, N. Miyazawa, K. Sharma, and P. K. Madhu, “Comparison among Magnus/Floquet/Fer expansion schemes in solid-state NMR.” J. Chem. Phys. 142, 134201 (2015).
  * Mananga and Charpentier (2016) E. S. Mananga and T. Charpentier, “On the Floquet-Magnus expansion: Applications in solid-state nuclear magnetic resonance and physics,” Phys. Rep. 609, 1–49 (2016).
  * Mananga and Charpentier (2018) E. S. Mananga and T. Charpentier, “Revisiting the applications of Floquet–Magnus and Fer expansion approaches in physics and solid-state NMR,” Int. J. Mod. Phys. B 32, 1850236–30 (2018).
  * Giscard and Bonhomme (2020) P.-L. Giscard and C. Bonhomme, “Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR,” Phys. Rev. Research 2, 023081 (2020).
  * Blanes  _et al._ (2009) S. Blanes, F. Casas, J. A. Oteo, and J. Ros, “The Magnus expansion and some of its applications,” Phys. Rep. 470, 151–238 (2009).
  * Mote, Agarwal, and Madhu (2016) K. R. Mote, V. Agarwal, and P. K. Madhu, “Five decades of homonuclear dipolar decoupling in solid-state NMR: Status and outlook,” Prog. NMR Spectr. 97, 1–39 (2016).
  * Ernst (2003) M. Ernst, “Heteronuclear spin decoupling in solid-state NMR under magic-angle sample spinning,” J. Magn. Reson. 162, 1–34 (2003).
  * Augustine (1995) M. P. Augustine,  _Theoretical treatment of time dependent phenomena and the investigation of highly magnetized systems - ProQuest_ , Ph.D. thesis, Ann Arbor (1995).
  * Boender (1996) G. J. Boender,  _The stacking of chlorophylls in chlorosomal antennae of green bacteria_ , Ph.D. thesis, Leiden University, Leiden (1996).
  * Boender, Vega, and de Groot (1998) G. Boender, S. Vega, and H. J. M. de Groot, “A physical interpretation of the Floquet description of magic angle spinning nuclear magnetic resonance spectroscopy,” Mol. Phys. 95, 921–934 (1998).
  * Ernst, Samoson, and Meier (2005) M. Ernst, A. Samoson, and B. H. Meier, “Decoupling and recoupling using continuous-wave irradiation in magic-angle-spinning solid-state NMR: A unified description using bimodal Floquet theory,” J. Chem. Phys. 123, 064102 (2005).
  * Ramachandran and Griffin (2005) R. Ramachandran and R. G. Griffin, “Multipole-multimode Floquet theory in nuclear magnetic resonance.” J. Chem. Phys. 122, 164502 (2005).
  * Scholz, Meier, and Ernst (2007) I. Scholz, B. H. Meier, and M. Ernst, “Operator-based triple-mode Floquet theory in solid-state NMR,” J. Chem. Phys. 127, 204504 (2007).
  * Van Vleck (1929) J. H. Van Vleck, “On $\sigma$-type doubling and electron spin in the spectra of diatomic molecules,” Phys. Rev. 33, 467–506 (1929).
  * Primas (1961) H. Primas, “Eine verallgemeinerte störungstheorie für quantenmechanische mehrteilchenprobleme,” Helv. Phys. Acta 34, 331–351 (1961).
  * Primas (1963) H. Primas, “Generalized Perturbation Theory in Operator Form,” Rev. Mod. Phys. 35, 710–& (1963).
  * Dyson (1949) F. J. Dyson, “The Radiation Theories of Tomonaga, Schwinger, and Feynman,” Phys. Rev. 75, 486–502 (1949).
  * Ramesh and Krishnan (2001) R. Ramesh and M. S. Krishnan, “Effective Hamiltonians in Floquet theory of magic angle spinning using van Vleck transformation,” J. Chem. Phys. 114, 5967–5973 (2001).
  * Vinogradov, Madhu, and Vega (2001) E. Vinogradov, P. K. Madhu, and S. Vega, “Phase modulated Lee-Goldburg magic angle spinning proton nuclear magnetic resonance experiments in the solid state: A bimodal Floquet theoretical treatment,” J. Chem. Phys. 115, 8983–9000 (2001).
  * Cornwell (1984) J. F. Cornwell,  _Group Theory in Physics_ (Academic Press, 1984).
  * Plancherel and Leffler (1910) M. Plancherel and M. Leffler, “Contribution à l’étude de la représentation d’une fonction arbitraire par des intégrales définies,” Rendiconti del Circolo Matematico di Palermo 30, 289–335 (1910).
  * Cartan (1966) E. Cartan, “The Theory of Spinors,” The MIT Press, Cambridge, MA (1966).
  * Blümich and Spiess (1985) B. Blümich and H. W. Spiess, “Quaternions as a Practical Tool for the Evaluation of Composite Rotations,” J. Magn. Reson. 61, 356–362 (1985).
  * Rodrigues (1840) O. Rodrigues, “Des lois géométriques qui régissent les déplacements d’un système solide dans l’espace, et de la variation des coordonnées provenant de ces déplacements …,” J. Math. 5, 380–440 (1840).
  * Scholz  _et al._ (2008) I. Scholz, M. Huber, T. Manolikas, B. H. Meier, and M. Ernst, “MIRROR recoupling and its application to spin diffusion under fast magic-angle spinning,” Chem. Phys. Lett. 460, 278–283 (2008).
  * Wittmann  _et al._ (2014) J. J. Wittmann, L. Hendriks, B. H. Meier, and M. Ernst, “Controlling spin diffusion by tailored rf-irradiation schemes,” Chem. Phys. Lett. 608, 60–67 (2014).
  * Wittmann  _et al._ (2016) J. J. Wittmann, V. Agarwal, J. Hellwagner, A. Lends, R. Cadalbert, B. H. Meier, and M. Ernst, “Accelerating proton spin diffusion in perdeuterated proteins at 100 kHz MAS,” J. Biomol. NMR. 66, 233–242 (2016).
  * Bartle (2014) R. G. Bartle,  _The Elements of Integration and Lebesgue Measure_ (John Wiley & Sons, 2014).
  * Veshtort and Griffin (2011) M. Veshtort and R. G. Griffin, “Proton-driven spin diffusion in rotating solids via reversible and irreversible quantum dynamics,” J. Chem. Phys. 135, 134509 (2011).
  * Smith  _et al._ (1994) S. A. Smith, T. O. Levante, B. H. Meier, and R. R. Ernst, “Computer-Simulations in Magnetic-Resonance - an Object-Oriented Programming Approach,” J. Magn. Reson. Ser. A 106, 75–105 (1994).
  * Samoson and Pines (1989) A. Samoson and A. Pines, “Double rotor for solid-state NMR,” Rev. Sci. Instrum. 60, 3239–3241 (1989).

##  V Supplementary Information

##  VI Derivation of Van Vleck perturbation theory on a continuous frequency
space

In the following we apply Van Vleck perturbation theory on the single-mode
Floquet Hamiltonian. As usual we split the Hamiltonian into two parts

$\displaystyle{\mathcal{H}}_{F}={\mathcal{H}}_{F}^{(0)}+\varepsilon\,{\mathcal{H}}_{F}^{(1)}$
(146)

Next we apply a unitary transformation called van Vleck transformation on Eq.
(146) which is defined as

$\displaystyle\Lambda_{F}=\textrm{U}\,{\mathcal{H}}_{F}\,\textrm{U}^{\dagger}={\mathcal{H}}_{F}^{(0)}+W_{F}$
(147)

where we choose $\textrm{U}=e^{S}$ such that
$[\Lambda_{F},{\mathcal{H}}_{F}^{(0)}]=0$. This commutation relation is
fulfilled if $[\textrm{U},{\mathcal{H}}_{F}^{(0)}]=0$. We define the nested
commutator as
$[S,{\mathcal{H}}_{F}^{(0)}]_{m}=[S,[S,{\mathcal{H}}_{F}^{(0)}]_{m-1}]$ with
$[S,{\mathcal{H}}_{F}^{(0)}]_{0}={\mathcal{H}}_{F}^{(0)}$. In this notation
the Baker-Campbell-Hausdorff formula takes the form

$\displaystyle\Lambda_{F}=\sum_{m=0}^{\infty}\frac{[S,{\mathcal{H}}_{F}]_{m}}{m!}={\mathcal{H}}^{(0)}_{F}+\sum_{n=1}^{\infty}\varepsilon^{n}\Lambda_{F}^{(n)}$
(148)

where we identify

$\displaystyle
W_{F}=\sum_{n=1}^{\infty}\varepsilon^{n}\Lambda_{F}^{(n)}=\sum_{m=1}^{\infty}\frac{[S,{\mathcal{H}}_{F}]_{m}}{m!}$
(149)

This is a important equation, because it connects the Hamiltonian
${\mathcal{H}}_{F}$ and the operator $S$ with the perturbation series. We now
proceed with inserting the expansion

$\displaystyle S=\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)}$ (150)

as well as Eq. (146) in Eq. (149), which leads to

$\displaystyle\sum_{n=1}^{\infty}\varepsilon^{n}\Lambda_{F}^{(n)}$
$\displaystyle=\sum_{m=1}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},{\mathcal{H}}_{F}^{(0)}+\varepsilon\,{\mathcal{H}}_{F}^{(1)}]_{m}}{m!}$
$\displaystyle=\sum_{m=1}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},{\mathcal{H}}_{F}^{(0)}]_{m}+[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},\varepsilon\,{\mathcal{H}}_{F}^{(1)}]_{m}}{m!}$
$\displaystyle=\sum_{m=1}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},{\mathcal{H}}_{F}^{(0)}]_{m}}{m!}+\sum_{m=1}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},\varepsilon\,{\mathcal{H}}_{F}^{(1)}]_{m}}{m!}$
$\displaystyle=\sum_{l=1}^{\infty}\varepsilon^{l}[S^{(l)},{\mathcal{H}}_{F}^{(0)}]+\sum_{m=2}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},{\mathcal{H}}_{F}^{(0)}]_{m}}{m!}+\sum_{m=1}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},\varepsilon\,{\mathcal{H}}_{F}^{(1)}]_{m}}{m!}$
(151)

As a result, we arrive at the equation

$\displaystyle[S^{(l)},{\mathcal{H}}_{F}^{(0)}]=\Lambda_{F}^{(l)}-\Phi^{(l)}_{F}$
(152)

with

$\displaystyle\sum_{j=1}^{\infty}\varepsilon^{j}\Phi^{(j)}_{F}=\sum_{m=2}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},{\mathcal{H}}_{F}^{(0)}]_{m}}{m!}+\sum_{m=1}^{\infty}\frac{[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},\varepsilon\,{\mathcal{H}}_{F}^{(1)}]_{m}}{m!}$
(153)

Utilizing mathematical induction over $m$ we can prove following identities

$\displaystyle\sum_{m=1}^{\infty}[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},\varepsilon{\mathcal{H}}_{F}^{(1)}]_{m}$
$\displaystyle=\sum_{m=1}^{\infty}\sum_{n_{1}=1}^{\infty}\sum_{n_{2}=1}^{\infty}\cdots\sum_{n_{m}=1}^{\infty}\,\varepsilon^{1+n_{1}+n_{2}+\cdots+n_{m}}[S^{(n_{1})},[S^{(n_{2})},\cdots[S^{(n_{m})},{\mathcal{H}}^{(1)}_{F}]\cdots]]$
(154)
$\displaystyle\sum_{m=2}^{\infty}[\sum_{l=1}^{\infty}\varepsilon^{l}S^{(l)},{\mathcal{H}}_{F}^{(0)}]_{m}$
$\displaystyle=\sum_{m=2}^{\infty}\sum_{n_{1}=1}^{\infty}\sum_{n_{2}=1}^{\infty}\sum_{n_{3}=1}^{\infty}\cdots\sum_{n_{m}=1}^{\infty}\varepsilon^{n_{1}+n_{2}+\cdots+n_{m}}[S^{(n_{1})},[S^{(n_{2})},\cdots[S^{(n_{m})},{\mathcal{H}}^{(0)}_{F}]\cdots]]$
(155)

which lead to a more explicit form of
$\sum_{j=1}^{\infty}\varepsilon^{j}\Phi^{(j)}_{F}$

$\displaystyle\sum_{j=1}^{\infty}\varepsilon^{j}\Phi^{(j)}_{F}$
$\displaystyle=\sum_{m=2}^{\infty}\sum_{n_{1}=1}^{\infty}\sum_{n_{2}=1}^{\infty}\cdots\sum_{n_{m}=1}^{\infty}\frac{\varepsilon^{n_{1}+n_{2}+\cdots+n_{m}}}{m!}[S^{(n_{1})},[S^{(n_{2})},\cdots[S^{(n_{m})},{\mathcal{H}}^{(0)}_{F}]\cdots]]$
$\displaystyle\hskip
14.22636pt+\sum_{m=1}^{\infty}\sum_{n_{1}=1}^{\infty}\sum_{n_{2}=1}^{\infty}\cdots\sum_{n_{m}=1}^{\infty}\,\frac{\varepsilon^{1+n_{1}+n_{2}+\cdots+n_{m}}}{m!}[S^{(n_{1})},[S^{(n_{2})},\cdots[S^{(n_{m})},{\mathcal{H}}^{(1)}_{F}]\cdots]]$
(156)

Although this expression is quite complicated, it allows relatively easy and
efficient calculation of $\Phi^{(l)}_{F}$. To find a solution of Eq. (152) for
$S^{(l)}$ we adapt the approach of Primas Primas (1961, 1963) and use
projection operator defined as

$\displaystyle\Pi(X)=\int P(u)XP(u)\,\textrm{d}u$
$\displaystyle{\mathcal{H}}^{(0)}_{F}=\int\limits_{-\infty}^{\infty}\kappa(u)P(u)\,\textrm{d}u$
(157)

By applying the projection operator $\Pi(X)$ on Eq. (152) we obtain

$\displaystyle\Pi([S^{(l)},{\mathcal{H}}_{F}^{(0)}])=\Pi(\Lambda_{F}^{(l)}-\Phi^{(l)}_{F})$
(158)

Since $[\Lambda_{F},{\mathcal{H}}_{F}^{(0)}]=0$, hence
$[\textrm{U},{\mathcal{H}}_{F}^{(0)}]=0$, we have

$\displaystyle
0=\Pi([S^{(l)},{\mathcal{H}}_{F}^{(0)}])=\Pi(\Lambda_{F}^{(l)}-\Phi^{(l)}_{F})=\Lambda_{F}^{(l)}-\Pi(\Phi^{(l)}_{F})\implies\Lambda_{F}^{(l)}=\Pi(\Phi^{(l)}_{F})$
(159)

By introducing the commutation operator

$\displaystyle\Gamma_{F}(X)=[X,{\mathcal{H}}_{F}^{(0)}],\hskip
28.45274pt\Gamma_{F}^{-1}(X)=\int\int\frac{P(u)XP(v)}{\kappa(u)-\kappa(v)}\,\textrm{d}u\,\textrm{d}v$
(160)

we can rewrite Eq. (152) as

$\displaystyle\Gamma_{F}(S^{(l)})=\Pi(\Phi^{(l)}_{F})-\Phi^{(l)}_{F}$ (161)

and finally arrive at the formal solution

$\displaystyle S^{(l)}=\Gamma_{F}^{-1}(\Pi(\Phi^{(l)}_{F})-\Phi^{(l)}_{F})$
(162)

The solution has the same form as found by Primas Primas (1961, 1963), but
with differ in the operators $\Pi(X)$ and $\Gamma_{F}(X)$ given in Eq.(157)
and Eq. (160), respectively. Ultimately we have all the ingredients to derive
each term of the perturbation series.

###  VI.1 Equivalence of the representations of the Floquet Hamiltonian (Eq.
(23) and Eq. (25) of the main text)

The following shows that Eq. (23) and Eq. (25) of the main text is equivalent:

$\displaystyle\bra{\chi,\mu}\mathcal{H}\ket{\xi,\nu}$
$\displaystyle=\bra{\chi,\mu}\int
D(\Omega)\otimes\widehat{{\mathcal{H}}}(\Omega)\,\textrm{d}\Omega+\hat{\Omega}\otimes\mathbf{1}\,\ket{\xi,\nu}$
$\displaystyle=\bra{\chi,\mu}\int
D(\Omega)\otimes\widehat{{\mathcal{H}}}(\Omega)\ket{\xi,\nu}\textrm{d}\Omega+\bra{\chi,\mu}\hat{\Omega}\otimes\mathbf{1}\,\ket{\xi,\nu}$
$\displaystyle=\int\bra{\chi}D(\Omega)\ket{\xi}\otimes\bra{\mu}\widehat{{\mathcal{H}}}(\Omega)\ket{\nu}\textrm{d}\Omega+\bra{\chi}\hat{\Omega}\ket{\xi}\otimes\bra{\mu}\mathbf{1}\ket{\nu}$
$\displaystyle=\int\braket{\chi}{\xi+\Omega}\otimes\bra{\mu}\widehat{{\mathcal{H}}}(\Omega)\ket{\nu}\textrm{d}\Omega+\xi\braket{\chi}{\xi}\otimes\braket{\mu}{\nu}$
$\displaystyle=\int\delta(\chi-\xi-\Omega)\otimes\widehat{{\mathcal{H}}}_{\mu\nu}(\Omega)\,\textrm{d}\Omega+\xi\,\delta(\chi-\xi)\otimes\delta_{\mu\nu}$
$\displaystyle=\mathbf{1}\otimes\widehat{{\mathcal{H}}}_{\mu\nu}(\chi-\xi)+\xi\,\delta(\chi-\xi)\otimes\delta_{\mu\nu}$
(163)

###  VI.2 Derivation of the first and second-order effective Hamiltonian for a
single-mode problem

Next we are going to use the previous results to calculate the first and
second-order effective Hamiltonian. For the sake of simplicity, we just
present the single mode case, because it already captures the whole procedure.
A generalization to the $n$-modal case in retrospect is straight forward,
therefore we will just state the results afterwards. The goal is to derive the
first three terms of the perturbation series, i.e.
$\bra{\Omega}{\mathcal{H}}^{(0)}_{F}+\Lambda^{(1)}_{F}+\Lambda^{(2)}_{F}\ket{\Omega^{\prime}}$.
In essence, the only equations we need are Eq. (156) (or Eq.(153)), Eq. (159)
and Eq. (162). For the single mode case the Hamiltonians are

$\displaystyle{\mathcal{H}}^{(0)}_{F}$ $\displaystyle=\widehat{\Omega}$ (164)
$\displaystyle{\mathcal{H}}^{(1)}_{F}$
$\displaystyle=\int\widehat{D}(\Omega)\otimes\widehat{{\mathcal{H}}}(\Omega)\,\textrm{d}\Omega=\Phi^{(1)}_{F}$
(165)

In the subsequent derivation the projection operator will be often applied on
the translation operator, which results in a delta distribution

$\displaystyle\Pi(\widehat{D}(\Omega))=\delta(\Omega)$ (166)

The projection of $\Phi^{(1)}_{F}$ simply leads to

$\displaystyle\Lambda_{F}^{(1)}=\Pi\left(\Phi^{(1)}_{F}\right)=\Pi\left(\int\widehat{D}(\Omega)\otimes\widehat{{\mathcal{H}}}(\Omega)\,\textrm{d}\Omega\right)=\int\delta(\Omega)\otimes\widehat{{\mathcal{H}}}(\Omega)\,\textrm{d}\Omega=\widehat{D}(0)\otimes\widehat{{\mathcal{H}}}(0)$
(167)

Inserting it into the Eq. (162) gives

$\displaystyle
S^{(1)}=\Gamma^{-1}\left(\Pi\left(\Phi^{(1)}_{F}\right)-\Phi^{(1)}_{F}\right)=\Gamma^{-1}\left(\widehat{D}(0)\otimes\widehat{{\mathcal{H}}}(0)-\int\widehat{D}(\Omega)\otimes\widehat{{\mathcal{H}}}(\Omega)\,\textrm{d}\Omega\right)=-\,PV\int\frac{\widehat{D}(\Omega)\otimes\widehat{{\mathcal{H}}}(\Omega)\,}{\Omega}\,\textrm{d}\Omega$
(168)

Notice that we use Cauchy principal value (PV) for the regularization of the
integral defined as

$\displaystyle PV\int_{a}^{c}f(x)=\lim\limits_{\epsilon\to
0^{+}}\left[\int_{a}^{b-\epsilon}f(x)\,dx+\int_{b+\epsilon}^{c}f(x)\,dx\right]$
(169)

Using this regularization technique is justified, since we subtract
$\widehat{D}(0)\otimes\widehat{{\mathcal{H}}}(0)$, which is the kernel of the
integral at the critical value $\Omega=0$ and therefore can be exclude it from
the integration. Equipped with the expression for $S^{(1)}$ we move to the
second order

$\displaystyle\Phi^{(2)}_{F}$
$\displaystyle=[S^{(1)},{\mathcal{H}}^{(1)}_{F}]+\frac{1}{2}[S^{(1)},[S^{(1)},{\mathcal{H}}_{F}^{(0)}]$
(170)
$\displaystyle=-\,PV\int\frac{[\widehat{{\mathcal{H}}}(\Omega),{\mathcal{H}}(0)]\widehat{D}(\Omega)}{\Omega}\,\textrm{d}\Omega+\frac{1}{2}\,PV\int\int\frac{[\widehat{{\mathcal{H}}}(\Omega),\widehat{{\mathcal{H}}}(\Omega^{\prime})]\widehat{D}(\Omega+\Omega^{\prime})}{\Omega}\,\textrm{d}\Omega\,\textrm{d}\Omega^{\prime}$
(171)

which results in

$\displaystyle\Lambda_{F}^{(2)}=\Pi\left(\Phi^{(2)}_{F}\right)=\frac{1}{2}\,PV\int\frac{[\widehat{{\mathcal{H}}}(\Omega),\widehat{{\mathcal{H}}}(-\Omega)]\widehat{D}(0)}{\Omega}\,\textrm{d}\Omega$
(172)

Finally we express the Hamiltonian in the product basis

$\displaystyle\bar{{\mathcal{H}}}^{(0)}+\bar{{\mathcal{H}}}^{(1)}+\bar{{\mathcal{H}}}^{(2)}=\bra{\Omega,\mu}{\mathcal{H}}_{F}^{(0)}+\Lambda^{(1)}_{F}+\Lambda^{(2)}_{F}\ket{\Omega^{\prime},\nu}$
(173)

with

$\displaystyle\bar{{\mathcal{H}}}^{(0)}=\widehat{\Omega}\hskip
28.45274pt\bar{{\mathcal{H}}}^{(1)}=\int\limits_{-\infty}^{\infty}\widehat{{\mathcal{H}}}(\Omega)\,\textrm{d}\Omega\hskip
28.45274pt\bar{{\mathcal{H}}}^{(2)}=\frac{1}{2}\,PV\int\int\frac{[\widehat{{\mathcal{H}}}(\Omega),\widehat{{\mathcal{H}}}(-\Omega)]}{\Omega}\textrm{d}\Omega$
(174)

As expected this results are similar to the corresponding equations known from
Floquet theory. The main difference is that we have an integral and a
translation operator instead of a sum and a ladder operator.
One of the main differences in the higher modal cases is the appearance of
resonance conditions between the modes. For the n-modal case the resonance
conditions are

$\displaystyle\Omega^{(0)}_{1}+\Omega^{(0)}_{2}+\cdots+\Omega^{(0)}_{n}=0$
(175)

The derivation of the effective Hamiltonians is quite similar and leads to the
first and second-order effective Hamiltonian, for the bimodal case

$\displaystyle\bar{{\mathcal{H}}}^{(1)}$
$\displaystyle=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\widehat{{\mathcal{H}}}(\Omega_{1},\Omega_{2})\,\textrm{d}\Omega_{1}\,\textrm{d}\Omega_{2}$
(176) $\displaystyle\bar{{\mathcal{H}}}^{(2)}$
$\displaystyle=\int\textrm{d}\Omega_{1}\int\textrm{d}\Omega_{1}^{(0)}\int\textrm{d}\Omega_{2}\int\textrm{d}\Omega_{2}^{(0)}\frac{[\widehat{{\mathcal{H}}}(\Omega_{1},\Omega_{2}),\widehat{{\mathcal{H}}}(\Omega_{1}^{(0)}-\Omega_{1},\Omega_{2}^{(0)}-\Omega_{2})]}{\Omega_{1}+\Omega_{2}}$
(177)

##  VII Calculation of the frequency-domain interaction-frame trajectory

###  VII.1 Propagator

The propagator of the $j$th time slice $\textrm{U}_{j}$ (Eq. (59)), has the
complex $2\times 2$ matrix representation

$\displaystyle\textrm{U}_{j}$
$\displaystyle=\exp\left(-i\,{\mathcal{H}}_{j}\delta
t\right)=\frac{1}{|\theta|\delta
t}\left(\begin{array}[]{cc}|\theta|\cos\left(\frac{|\theta|\delta
t}{2}\right)+i\,\theta_{z}\sin\left(\frac{|\theta|\delta
t}{2}\right)&\sin\left(\frac{|\theta|\delta
t}{2}\right)(\theta_{y}+i\,\theta_{x})\\\ \sin\left(\frac{|\theta|\delta
t}{2}\right)(-\theta_{y}+i\,\theta_{x})&|\theta|\cos\left(\frac{|\theta|\delta
t}{2}\right)-i\,\theta_{z}\sin\left(\frac{|\theta|\delta t}{2}\right)\\\
\end{array}\right)$ (180)

where we dropped the index $j$ of $\theta$ for sake of simplicity. This
representation is used for derivation of Eq. (61) in the main text.

###  VII.2 The frequency-domain interaction-frame trajectory for a cyclic
pulse scheme

The interaction-frame trajectory of a cyclic pulse scheme with an basic
element of duration $\tau$ and an overall duration $N\tau$ can be written as

$\displaystyle a_{\mu\nu}(t)=\begin{dcases}a_{\mu\nu}(t)&0\leq t\leq\tau\\\
\sum_{\chi}a_{\mu\chi}(t-\tau)a_{\chi\nu}(\tau)&\tau\leq t\leq 2\tau\\\
\sum_{\chi}a_{\mu\chi}(t-2\tau)a_{\chi\nu}(2\tau)&2\tau\leq t\leq 3\tau\\\
\hskip 56.9055pt\vdots&\quad\quad\vdots\\\
\sum_{\chi}a_{\mu\chi}(t-(N-1)\tau)a_{\chi\nu}((N-1)\tau)&(N-1)\tau\leq t\leq
N\tau\end{dcases}$ (181)

Fourier transformation leads to

$\displaystyle\widehat{a}_{\mu\nu}(\Omega)$ $\displaystyle=\frac{1}{2\pi
N\tau}\biggl{(}\sum_{\chi}\int_{0}^{\tau}a_{\mu\nu}(t)e^{-i\Omega
t}dt+\int_{\tau}^{2\tau}a_{\mu\chi}(t-\tau)a_{\chi\nu}(\tau)e^{-i\Omega
t}dt+\int_{2\tau}^{3\tau}a_{\mu\chi}(t-2\tau)a_{\chi\nu}(2\tau)e^{-i\Omega
t}dt$ (182)
$\displaystyle\quad\quad+\cdots+\int_{(N-1)\tau}^{N\tau}a_{\mu\chi}(t-(N-1)\tau)a_{\chi\nu}((N-1)\tau)e^{-i\Omega
t}dt\biggr{)}$ $\displaystyle=\frac{1}{2\pi
N\tau}\sum_{\chi}\sum_{n=1}^{N}\int_{(n-1)\tau}^{n\tau}a_{\mu\chi}(t-(n-1)\tau)a_{\chi\nu}((n-1)\tau)e^{-i\Omega
t}dt$ (183) $\displaystyle=\frac{1}{2\pi
N\tau}\sum_{\chi}\sum_{n=1}^{N}e^{-i\Omega(n-1)\tau}\int_{0}^{\tau}a_{\mu\chi}(t^{\prime})a_{\chi\nu}((n-1)\tau)e^{-i\Omega
t^{\prime}}dt^{\prime}$ (184) $\displaystyle=\frac{1}{2\pi
N\tau}\sum_{\chi}\int_{0}^{\tau}a_{\mu\chi}(t^{\prime})e^{-i\Omega
t^{\prime}}dt^{\prime}\sum_{n=0}^{N-1}e^{-i\Omega n\tau}a_{\chi\nu}(n\tau)$
(185)
$\displaystyle=\sum_{\chi}\bar{a}_{\mu\chi}(\Omega)\,g^{(N)}_{\chi\nu}(\Omega)$
(186)

and therefore

$\displaystyle
g^{(N)}_{\mu\nu}(\Omega)\coloneqq\frac{1}{N}\sum_{n=0}^{N-1}e^{-in\Omega\tau}a_{\mu\nu}(n\tau)$
(187)

Notice that $a_{\mu\nu}(n\tau)=[\mathbf{a}^{n}(\tau)]_{\mu\nu}$ where
$\mathbf{a}(\tau)$ is the matrix with the elements $a_{\mu\nu}(\tau)$.
Furthermore we have

$\displaystyle\bar{a}_{\mu\nu}(\Omega)$
$\displaystyle=\frac{1}{2\pi\tau}\int\limits_{0}^{\tau}{a}_{\mu\nu}(t)e^{-i\Omega
t}dt$ (188)
$\displaystyle=\mathcal{F}\\{\widetilde{a}_{\mu\nu}(t)\Pi\left(t/\tau-1/2\right)\\}$
(189)
$\displaystyle=\mathcal{F}\\{\widetilde{a}_{\mu\nu}(t)\\}*\mathcal{F}\left\\{\Pi(t/\tau-1/2)\right\\}$
(190)

and

$\displaystyle\mathcal{F}\\{\widetilde{a}_{\mu\nu}(t)\\}=\widetilde{a}_{\mu\nu}(\Omega)=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}{a}_{\mu\nu}(t\,\textrm{mod}\tau)e^{-i\Omega
t}dt$ (191)

Notice that $\widetilde{a}_{\mu\nu}(\Omega)$ is just calculated using FFT. As
described in the main text, the special case of a cyclic interaction-frame
trajectory occurs when $\mathbf{a}(\tau)=\mathbf{1}$.