Visualization of dynamics in coupled multi-spin systems

Since the dawn of quantum mechanics, ways to visualize spins and their interactions have attracted the attention of researchers and philosophers of science. In this work we present a generalized measurement-based 3D-visualization approach for describing dynamics in strongly coupled spin ensembles. The approach brings together angular momentum probability surfaces (AMPS), Husimi Q functions, and DROPS (discrete representations of operators for spin systems) and finds particular utility when the total angular momentum basis is used for describing Hamiltonians. We show that, depending on the choice of a generalized measurement operator, the plotted surfaces either represent probabilities of finding the maximal projection of an angular momentum along any direction in space or represent measurable coherences between the states with different total angular momenta. Such effects are difficult to grasp by looking at (time-dependent) numerical values of density-matrix elements. The approach is complete in a sense that there is one-to-one correspondence between the plotted surfaces and the density matrix. Three examples of nuclear spin dynamics in two-spin systems are visualized: (i) a zero- to ultralow-field (ZULF) nuclear magnetic resonance (NMR) experiment in the presence of a magnetic field applied perpendicularly to the sensitive axis of the detector, (ii) interplay between chemical exchange and spin dynamics during high-field signal amplification by reversible exchange (SABRE), and (iii) a high-field spin-lock-induced crossing (SLIC) sequence, with the initial state being the singlet state between two spins. The presented visualization technique facilitates intuitive understanding of spin dynamics during complex experiments as exemplified here by the considered cases. Temporal sequences (“the movies”) of such surfaces show phenomena like interconversion of spin order between the coupled spins and are particularly relevant in ZULF NMR.

1 Introduction

The evolution of spins in nuclear magnetic resonance (NMR) experiments can be highly complex and non-intuitive. Thus, approaches to visualize spin dynamics of nuclear multi-spin systems are sought for both communication and research purposes. Most NMR textbooks visualize the motion of a single spin-1/2 (or an ensemble of spins) using the Bloch vector (e.g., Bloch, 1946; Feynman et al., 1957; Levitt, 2013).

In more complex cases such as coupled spin systems, visualization is significantly less straightforward, and the discussion is often assisted by drawing energy-level diagrams. These diagrams do not provide dynamic information but can be useful for representing populations of various spin states (e.g., Sørensen et al., 1984; Messiah, 1999; Levitt, 2013; Barskiy et al., 2019). Spin dynamics can be visually presented via coherence-transfer pathways, i.e., graphical conventions that show the path through different levels of coherences as a function of time during the implementation of high-field NMR pulse sequences. These pictures can be used to derive phase cycles and to obtain pure-phase NMR signals via symmetry properties of the phases appreciated from the graphical representations (Bodenhausen et al., 1984).

Recently, Garon et al. (2015) introduced an approach of so-called “DROPS” (DROPS: discrete representation of operators for spin systems) for visualizing spin operators (Garon et al., 2015; Leiner et al., 2017; Leiner and Glaser, 2018). This visualization approach is based on plotting 3D colored shapes (droplets) that correspond to linear combinations of spherical harmonics. The DROPS approach can represent interactions in Hamiltonians and in propagators as well as states of the density matrix (Garon et al., 2015; Leiner et al., 2017; Leiner and Glaser, 2018). The DROPS approach with a LISA basis (with defined linearity, subsystem, and auxiliary criteria) reflects the rotational symmetry of individual spins and can be useful for understanding high-field NMR experiments where each spin can be addressed individually (for example, when a strong magnetic field breaks the isotropic symmetry of the system). While the DROPS approach could be generalized to isotropic systems by using multiple tensor operator basis, it is challenging to extract the information corresponding to measurable properties from the DROPS with complicated colors. To address these limitations, we introduce a generalized approach based on angular momentum probability surfaces (AMPS) (Auzinsh, 1997; Rochester and Budker, 2001; Auzinsh et al., 2010).

The AMPS were introduced to visualize the angular momentum state of atoms. It is worth pointing out that AMPS are a particular case of the s-parametrized phase-space functions with $s=-\mathrm{1}$, also known as the Husimi Q function first published in 1940; see Appendix G (Husimi, 1940; Stratonovich, 1957; Agarwal, 1981; Várilly and Gracia-Bondía, 1989; Brif and Mann, 1999; Koczor et al., 2020).

To understand the physics behind the plotting, one can assume that the spin ensemble is measured by a Stern–Gerlach experiment, and the probability of finding the maximal projection state is plotted as a surface; the radius from the origin of the coordinate system to a point on the surface represents the measured maximal probability along that direction. The “measurements” are conducted with the measurement device rotated by Euler angles with the z–y–z convention $(\mathit{\varphi },\mathit{\theta },\mathrm{0})$. The AMPS can be constructed by setting the radius of the surface along the direction (θ,ϕ) equal to the measured probability. Despite their utility in atomic physics, AMPS have limitations for use in NMR. For example, they do not represent coherences connecting states with different total angular momenta (including states of the same total angular momentum quantum number F but belonging to different manifolds) which play a crucial role in NMR of multi-spin systems. In this study, we demonstrate that such coherences can also be visualized assuming a generalized form of the measurement operator, with the measurement device rotated along various directions. In addition, we show that one of the selected measurements is directly related to the measured zero- to ultralow-field (ZULF) NMR signal. As a result, our approach is helpful in understanding the complicated splittings in ZULF NMR spectra. The approach presented in this work constitutes convenient means for visualizing complex dynamics in multi-spin systems exemplified here for pairs of nuclear spins-1/2.

2 Results

In our visualization approach, a density matrix $\widehat{\mathit{\rho }}$ is represented by the action of a measurement operator ${\widehat{O}}_{\widehat{n}}$:

$$\begin{array}{}\text{(1)}& r_{\widehat{n}}=\mathrm{Tr}\left(\widehat{\mathit{\rho }}{\widehat{O}}_{\widehat{n}}\right).\end{array}$$

Here, $r_{\widehat{n}}$ indicates the distance of the plotted surface point to the origin along $\widehat{n}$ (Fig. 1a). In order to find the exact form of the operator ${\widehat{O}}_{\widehat{n}}$, a measurement is first defined as ${\widehat{O}}_{\widehat{z}}$ (note that the subscript $\widehat{z}$ denotes the direction of plotting the result of the measurement and does not specify the direction of the measured property). For example, ${\widehat{O}}_{\widehat{z}}$ could be $({\widehat{I}}_{\mathrm{1}x}{\widehat{I}}_{\mathrm{2}y}-{\widehat{I}}_{\mathrm{1}y}{\widehat{I}}_{\mathrm{2}x})$; see below. The operator ${\widehat{O}}_{\widehat{n}}$ then represents the observable obtained via global rotation ${\widehat{R}}_{\widehat{z}\to \widehat{n}}$, which brings $\widehat{z}$ to $\widehat{n}$:

$$\begin{array}{}\text{(2)}& {\widehat{O}}_{\widehat{n}}={\widehat{R}}_{\widehat{z}\to \widehat{n}}{\widehat{O}}_{\widehat{z}}{\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}}.\end{array}$$

Figure 1 (a) The visualized surface is a collection of points plotted at a distance $r_{\widehat{n}}$ from the origin along the direction $\widehat{n}$ (red color represents a positive number and blue color represents a negative number). (b) Density matrix written on a total angular momentum basis. In such a basis, the density operator is decomposed into blocks according to F (total angular momentum quantum number); the diagonal block is denoted as (F,F) and the off-diagonal block is denoted as (F,K). In the following figures we keep the same axis angles.

There is one important caveat. The observable ${\widehat{O}}_{\widehat{z}}$ must remain unchanged under the rotation about $\widehat{z}$; otherwise, the way in which the rotation is performed will affect the plotted result. To avoid ambiguity, only zero-quantum operators can be used as the measurement operators since they fulfill the required property of invariance under rotation (see Appendix A). Moreover, the operator ${\widehat{O}}_{\widehat{n}}$ must be Hermitian for a measurement $r_{\widehat{n}}$ to be real. Under such constraints, one may still obtain a negative value for $r_{\widehat{n}}$; thus, color is introduced as an additional “degree of freedom”. As a convention, in this work we set the distance of the surface along $\widehat{n}$ as the amplitude of $r_{\widehat{n}}$ and set the color of the surface along $\widehat{n}$ to be red if $r_{\widehat{n}}$ is positive and blue if $r_{\widehat{n}}$ is negative. If the observable operators are not set to be Hermitian, a color map is necessary to express (0, 2π) phase dependence of the calculated expectation value, which is, generally, a complex number (Garon et al., 2015).

In the following discussion, we focus on constructing measurement operators $\mathit{\left\{}{\widehat{O}}_{\widehat{n}}\mathit{\right\}}$ for isotropic or nearly isotropic coupled nuclear spin systems. Nuclear spin Hamiltonians for such systems remain unchanged (or changed insignificantly) under global rotation. Therefore, we choose to work in the total angular momentum basis which is defined by the eigenstates of the operators ${\widehat{F}}_z$ and $(\widehat{\mathbf{F}},\widehat{\mathbf{F}})={\widehat{F}}_x^{\mathrm{2}}+{\widehat{F}}_y^{\mathrm{2}}+{\widehat{F}}_z^{\mathrm{2}}={\widehat{F}}^{\mathrm{2}}$:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\widehat{F}}_z|F,m\rangle =m|F,m\rangle ,\\ & {\widehat{F}}^{\mathrm{2}}|F,m\rangle =F(F+\mathrm{1}\left)\right|F,m\rangle ,\end{array}\end{array}$$

where ${\widehat{F}}_{\mathit{\alpha }}$ is a total spin projection operator ($\mathit{\alpha }=x,y,z$), ${\widehat{F}}_{\mathit{\alpha }}={\sum }_{i=\mathrm{1}}^N{\widehat{I}}_{i\mathit{\alpha }}$ for N spins-1/2, and $\widehat{\mathbf{F}}=({\widehat{F}}_x,{\widehat{F}}_y,{\widehat{F}}_z)$ is the total angular momentum vector operator.

Now we write the density operator in the total angular momentum basis which decomposes it into blocks according to the total angular momentum quantum number F (Fig. 1b). Since blocks of the density matrix transform differently under global rotations, visualizations (Eq. 1) of different blocks are introduced separately to capture full information about the system. We denote diagonal blocks as (F,F) and non-diagonal blocks as (F,K). For the system consisting of N nuclear spins, there are in total $h=C_N^{N/\mathrm{2}}$ (binomial coefficient) diagonal blocks for the even N and $h=C_N^{(N-\mathrm{1})/\mathrm{2}}$ diagonal blocks for the odd N. In order to reproduce the full spin dynamics of a general density matrix using our method, a total number of h² visualized surfaces is required. In practice, depending on the specific initial state and interactions in the system, a lower number of visualizations is sufficient to fully represent dynamics since some blocks of the density matrix remain unpopulated during the experiment.

Since we consider Hermitian operators that are invariant under rotations about $\widehat{z}$ (termed hereafter zero-quantum operators and denoted ${\widehat{\mathrm{ZQ}}}_{\mathit{\phi },m}^{(F,K)}$), a generalized measurement operator in a laboratory frame is written as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\widehat{O}}_{\widehat{z}}& ={\widehat{\mathrm{ZQ}}}_{\mathit{\phi },m}^{(F,K)}\\ & ={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\left(e^{i\mathit{\phi }}|F,m\rangle \langle K,m|+e^{-i\mathit{\phi }}|K,m\rangle \langle F,m|\right),\end{array}\end{array}$$

where the operators are defined up to a phase φ and a projection quantum number m such that $\left|m\right|\le min(F,K)$.

When it comes to a diagonal block (F,F), one can choose a measurement operator ${\widehat{\mathrm{ZQ}}}_{\mathit{\pi }/\mathrm{4},F}^{(F,F)}$, which is the same as the maximum-projection operator $|F,F\rangle \langle F,F|$. The function of the surface representing the (F,F) block can thus be expressed using Eq. (1):

$$\begin{array}{}\text{(5)}& {}_{\mathit{\pi }/\mathrm{4},F}^{(F,F)}{r}_{\widehat{n}}=\mathrm{Tr}\left[\widehat{\mathit{\rho }}{\widehat{R}}_{\widehat{z}\to \widehat{n}}|F,F\rangle \langle F,F|{\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}}\right].\end{array}$$

One may see that up to this point the presented visualization method is equivalent to the previously described AMPS approach (Rochester and Budker, 2001; Auzinsh et al., 2010); ${}_{\mathit{\pi }/\mathrm{4},F}^{(F,F)}{r}_{\widehat{n}}$ presented in such a way is also known in the literature as the Husimi Q function (see Appendix G). Note that the visualization is complete in a sense that it fully represents a density matrix for block (F,F) (see Appendix D).

For the representation of a pair of off-diagonal blocks (F,K) and (K,F), in Eq. (4) we choose two measurement operators ${\widehat{\mathrm{ZQ}}}_{\mathit{\phi },m}^{(F,K)}$ such that φ=0 and $\mathit{\phi }=\mathit{\pi }/\mathrm{2}$ (with m chosen from the allowed range). Indeed, given a fixed m in Eq. (4), one may notice that the defined operators vary only by one number, φ, thus forming a 2D real operator space. Therefore, it is enough to consider only two values of φ to fully represent the spin dynamics of the coherences in blocks (F,K) and (K,F).

One can write the function of the surfaces according to Eq. (1):

$$\begin{array}{}\text{(6)}& {}_{\mathit{\phi },m}^{(F,K)}{r}_{\widehat{n}}=\mathrm{Tr}\left[\widehat{\mathit{\rho }}{\widehat{R}}_{\widehat{z}\to \widehat{n}}{\widehat{\mathrm{ZQ}}}_{\mathit{\phi },m}^{(F,K)}{\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}}\right],\end{array}$$

with φ being equal to 0 or $\mathit{\pi }/\mathrm{2}$. Similarly to AMPS, the two off-diagonal blocks (F,K) and (K,F) can be fully represented by these two surfaces (see Appendix D). To summarize, the proposed visualization approach is complete in a sense that there is a one-to-one correspondence between the collection of surfaces and blocks of the density operator.

As an example, consider a pair of spin-1/2 nuclei. The density matrix written in the total angular momentum basis consists of four blocks: (1,1), (1,0), (0,1), and (0,0). For the two diagonal blocks (1,1) and (0,0), one can consider measurements $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ and $\left(\right|\mathrm{0},\mathrm{0}\rangle \langle \mathrm{0},\mathrm{0}|{)}_{\widehat{n}}$. The two remaining off-diagonal blocks can be represented using the operators ${\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}$ and ${\left({\widehat{\mathrm{ZQ}}}_{\frac{\mathit{\pi }}{\mathrm{2}},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$. Only m=0 is permissible due to the requirement $\left|m\right|\le min(\mathrm{1},\mathrm{0})$. Thereby, one can construct the first measurement operator ${\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}=\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}({\widehat{I}}_{\mathrm{1}z}-{\widehat{I}}_{\mathrm{2}z})$ (i.e., the magnetization difference or the so-called zero-quantum in-phase coherence) when φ equals zero and the second operator ${\widehat{\mathrm{ZQ}}}_{\frac{\mathit{\pi }}{\mathrm{2}},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}=\sqrt{\mathrm{2}}({\widehat{I}}_{\mathrm{1}y}{\widehat{I}}_{\mathrm{2}x}-{\widehat{I}}_{\mathrm{1}x}{\widehat{I}}_{\mathrm{2}y})$ when φ equals $\mathit{\pi }/\mathrm{2}$ (the so-called zero-quantum coherence out of phase; Pravdivtsev et al., 2017).

Examples of different visualized surfaces are presented in Fig. 2. The first example shown in Fig. 2a is a polarized state that exemplifies magnetic orientation, i.e., a preferred direction of magnetization in space (see Appendix F for a step-by-step plotting tutorial). The second example in Fig. 2b shows an anti-phase spin order exemplifying alignment (Auzinsh et al., 2010). Such a spin order is often an initial state in parahydrogen-based hyperpolarization experiments at a high field (Bowers and Weitekamp, 1987). These density operators only have nonzero entries within the (1,1) and (0,0) blocks; therefore, the surfaces based on $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ measurements fully represent the related density operators (trivially, the (0,0) block is represented by a sphere with a radius $\mathrm{1}/\mathrm{4}$ and therefore is not shown). In Fig. 2c–d, we give an example of a density operator $\widehat{\mathit{\rho }}=(\mathrm{1}/\mathrm{4})\text{𝟙}+(\mathrm{1}/\mathrm{4})({\widehat{I}}_{\mathrm{1}z}-{\widehat{I}}_{\mathrm{2}z})$ that is often encountered in various experiments involving two-spin systems (see the discussion below). Figure 2c shows the visualized surface based on the measurement $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$, which appears isotropic. However, the state is intrinsically anisotropic along the $\widehat{z}$ axis, which can be displayed by using a measurement operator ${\left({\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$ (Fig. 2d). Through these examples, the usefulness of our visualization approach becomes apparent: one can spot the presence of different symmetries (or lack of thereof) by simply looking at the calculated surface without analyzing the explicit structure of the density matrices.

Figure 2 (a–c) Visualizations using the measurement observable $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ for the following density operators: (a) $\widehat{\mathit{\rho }}=(\mathrm{1}/\mathrm{4})\text{𝟙}+(\mathrm{1}/\mathrm{4})({\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z})$, representing a state orientated along $\widehat{z}$; (b) $\widehat{\mathit{\rho }}=(\mathrm{1}/\mathrm{4})\text{𝟙}+{\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}z}$, representing a state aligned along $\widehat{z}$; (c) $\widehat{\mathit{\rho }}=(\mathrm{1}/\mathrm{4})\text{𝟙}+(\mathrm{1}/\mathrm{4})({\widehat{I}}_{\mathrm{1}z}-{\widehat{I}}_{\mathrm{2}z})$, representing a state that appears isotropic from the point of measurement of the maximal projection of the total angular momentum with F=1. This state is intrinsically anisotropic, which is shown in (d) by plotting the surface with the measurement operator ${\left({\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$.

3 Discussion

3.1 Equivalence and symmetries

As discussed before, one can easily spot the symmetries of the density operator by looking at corresponding visualized surfaces (Fig. 2). Furthermore, as shown in Appendix D, the functions $r_{\widehat{n}}$ corresponding to surfaces constructed via Eqs. (5) and (6) contain as much information as the density matrix.

Figure 3 The density operators written in Zeeman basis {|αα〉, |αβ〉, |βα〉 and |ββ〉} (top) and their visualizations (bottom) for (a) $\widehat{\mathit{\rho }}=\frac{\mathrm{1}}{\mathrm{4}}\text{𝟙}+{\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}z}$, (b) the same density operator as in panel (a) rotated by $-\frac{\mathit{\pi }}{\mathrm{4}}$ around $\widehat{x}$, and (c) $\widehat{\mathit{\rho }}=\frac{\mathrm{1}}{\mathrm{4}}\text{𝟙}+\frac{\mathrm{1}}{\mathrm{4}}({\widehat{I}}_{\mathrm{1}z}-{\widehat{I}}_{\mathrm{2}z})$, (d) the same density operator as in panel (c) rotated by $-\frac{\mathit{\pi }}{\mathrm{4}}$ around $\widehat{x}$. Note that the surfaces in panels (a)–(b) and panels (c)–(d) are plotted using the measurement operators $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ and ${\left({\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$, respectively. The rotation is apparent by looking at the plotted surface, while it is hard to understand directly from the matrix form of the density operator.

Now we show explicitly that any global rotation applied to nuclear spins is directly reflected by the rotations of the plotted surfaces. Examples are given in Fig. 3. A state with alignment (see Fig. 2b), $\widehat{\mathit{\rho }}=(\mathrm{1}/\mathrm{4})\text{𝟙}+{\widehat{I}}_{\mathrm{1}z}{\widehat{I}}_{\mathrm{2}z}$, is shown in Fig. 3a; the density operator is fully represented by the surface when plotted using the measurement observable $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$. Figure 3b shows the same density rotated by $-(\mathit{\pi }/\mathrm{4})$ along $\widehat{x}$. As expected, the visualized surface clearly represents a rotated aligned state. The same logic also applies to the state $\widehat{\mathit{\rho }}=(\mathrm{1}/\mathrm{4})\text{𝟙}+(\mathrm{1}/\mathrm{4})({\widehat{I}}_{\mathrm{1}z}-{\widehat{I}}_{\mathrm{2}z})$ (Fig. 3c–d), for which the visualized surface plotted with $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ appears isotropic. Since the state is clearly anisotropic, another measurement operator should be chosen for visualization. A proper choice is ${\left({\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$; the symmetry of the surface is clearly seen when rotation is applied to the density operator (Fig. 3c–d).

Figure 4 Examples illustrating the relationship between the symmetry of the visualized surface and the elements of the density operator. (a) Density operator written in Zeeman basis plotted with $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$; (b) density operator written in Zeeman basis plotted with (left) $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ and (right) ${\left({\widehat{\mathrm{ZQ}}}_{\frac{\mathit{\pi }}{\mathrm{2}},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$. In panel (a) the surface is 2-fold symmetric around $\widehat{z}$, and the corresponding density operator only has coherences with $\left|\mathrm{\Delta }m\right|=\mathrm{0}$ and $\left|\mathrm{\Delta }m\right|=\mathrm{2}$. In panel (b), from left to right, the two surfaces are 2-fold symmetric and rotational symmetric (q-fold symmetric for an arbitrary positive integer q) around $\widehat{z}$, and the related density operator (written in Zeeman basis) only has coherences with $\left|\mathrm{\Delta }m\right|=\mathrm{0}$ and $\left|\mathrm{\Delta }m\right|=\mathrm{2}$.

In addition, the symmetry of the constructed surfaces can reflect useful information about elements of the density matrix. If a q-fold symmetry is present for all the surfaces, the density matrix, when written with the quantization axis set parallel to the symmetry axis, only has nonzero elements with Δm=qN, where N is an integer (see Appendix E). Examples are given in Fig. 4. In Fig. 4a, the only surface-plotted $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ representing the state is 2-fold symmetric along $\widehat{z}$. The related density matrix (written in Zeeman basis) only has nonzero coherences with $\left|\mathrm{\Delta }m\right|=\mathrm{0}$ and $\left|\mathrm{\Delta }m\right|=\mathrm{2}$. Figure 4b shows an example where two surfaces are needed to completely represent the state. The one plotted with the measurement operator $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ (left) is 2-fold symmetric around $\widehat{z}$, while that plotted with ${\left({\widehat{\mathrm{ZQ}}}_{\frac{\mathit{\pi }}{\mathrm{2}},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$ (middle) is rotationally symmetric (q-fold symmetric for an arbitrary positive integer q) around $\widehat{z}$. As a result, both surfaces possess a 2-fold symmetry around $\widehat{z}$, and the related density matrix only has nonzero coherences with $\left|\mathrm{\Delta }m\right|=\mathrm{0}$ and $\left|\mathrm{\Delta }m\right|=\mathrm{2}$.

3.2 Zero- to ultralow-field nuclear magnetic resonance

First, we show how one could use the presented visualization approach to better understand spin dynamics in ZULF NMR experiments (Ledbetter et al., 2011; Blanchard et al., 2021). As an example, consider the ¹³C-labeled formic acid with ¹H and ¹³C spins initially polarized along the z axis (direction of a magnetometer's sensitive axis). During the acquisition a weak bias magnetic field, B_x (perpendicular to the magnetometer's sensitive axis), is applied (Fig. 5a). In such an experiment, nuclear spins evolve due to the heteronuclear J coupling and the bias field, and the ZULF NMR spectrum is collected by measuring the magnetization along the z direction. A simulated NMR spectrum following the described procedure is shown in Fig. 5b. First, a low-frequency peak is positioned at the average Larmor frequency between the proton and the carbon, $\stackrel{\mathrm{‾}}{\mathit{\nu }}=({\mathit{\nu }}_{{}^{\mathrm{1}}\mathrm{H}}+{\mathit{\nu }}_{{}^{\mathrm{13}}\mathrm{C}})/\mathrm{2}$ (${\mathit{\nu }}_{{}^{\mathrm{1}}\mathrm{H}}$ and ${\mathit{\nu }}_{{}^{\mathrm{13}}\mathrm{C}}$ are the ¹H and ¹³C nuclear Larmor frequencies, respectively). In addition, a doublet with splitting equal to the sum of the Larmor frequencies of proton and carbon is centered at the J-coupling frequency.

Figure 5 Visualizations of spin dynamics in an AX system (¹H–¹³C nuclear pair) during the zero- to ultralow-field nuclear magnetic resonance (ZULF NMR) experiment. (a) Scheme of an exemplary ZULF NMR experiment in which a perpendicular field of 2.64 mG is applied during the detection. (b) Simulation of the corresponding ZULF NMR spectrum. Assume the initial density operator of the system is ${\widehat{\mathit{\rho }}}_{\mathrm{0}}=(\mathrm{1}/\mathrm{4})\text{𝟙}+(\mathrm{1}/\mathrm{10})({\widehat{I}}_{\mathrm{1}z}+\mathrm{4}{\widehat{I}}_{\mathrm{2}z})$ (${\widehat{\mathbf{I}}}_{\mathrm{1}}$ and ${\widehat{\mathbf{I}}}_{\mathrm{2}}$ denote ¹³C and ¹H nuclei, respectively). Panels (c)–(d) show surfaces representing spin evolution in the ZULF experiment plotted with the measurements (c) $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ over a timescale ${\mathit{\tau }}_{\mathrm{1}}=\mathrm{1}/\stackrel{\mathrm{‾}}{\mathit{\nu }}=\mathrm{2}/({\mathit{\nu }}_{{}^{\mathrm{1}}\mathrm{H}}+{\mathit{\nu }}_{{}^{\mathrm{13}}\mathrm{C}})$. Panel (d) ${\left({\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$ over a timescale τ₁ (inset) shows the evolution of the surfaces plotted in panel (d) over a shorter timescale ${\mathit{\tau }}_{\mathrm{2}}=\mathrm{1}/J$.

The observable operator describing a ZULF measurement is proportional to ${\mathit{\gamma }}_{{}^{\mathrm{13}}\mathrm{C}}{\widehat{I}}_{\mathrm{1}z}+{\mathit{\gamma }}_{{}^{\mathrm{1}}\mathrm{H}}{\widehat{I}}_{\mathrm{2}z}$ (the operators ${\widehat{\mathbf{I}}}_{\mathrm{1}}$ and ${\widehat{\mathbf{I}}}_{\mathrm{2}}$ denote the ¹³C and ¹H spins, respectively, and their projections are accordingly ${\widehat{I}}_{\mathrm{1}z}$ and ${\widehat{I}}_{\mathrm{2}z}$). Further, one can decompose the observable operator into a symmetric part $\left(({\mathit{\gamma }}_{{}^{\mathrm{13}}\mathrm{C}}+{\mathit{\gamma }}_{{}^{\mathrm{1}}\mathrm{H}})/\mathrm{2}\right)({\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z})$ and an antisymmetric part $\left(({\mathit{\gamma }}_{{}^{\mathrm{13}}\mathrm{C}}-{\mathit{\gamma }}_{{}^{\mathrm{1}}\mathrm{H}})/\mathrm{2}\right)({\widehat{I}}_{\mathrm{1}z}-{\widehat{I}}_{\mathrm{2}z})$. One may notice that symmetric and antisymmetric parts are proportional to the rank-1 component of the measurement operator $|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|$ and the measurement operator ${\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}$, respectively. As a result, the ZULF spectra could be understood by checking the evolution of the surfaces plotted with $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$ and ${\left({\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$ measurement operators separately as is done below.

Let us first examine the symmetric part of the measurement operator in the ZULF NMR experiment. Figure 5c visualizes the evolution of the density operator ${\widehat{\mathit{\rho }}}_{\mathrm{0}}=(\mathrm{1}/\mathrm{4})\text{𝟙}+(\mathrm{1}/\mathrm{10})({\widehat{I}}_{\mathrm{1}z}+\mathrm{4}{\widehat{I}}_{\mathrm{2}z})$ (we assume high polarizations and set them proportional to gyromagnetic ratios using ${\mathit{\gamma }}_{{}^{\mathrm{1}}\mathrm{H}}/{\mathit{\gamma }}_{{}^{\mathrm{13}}\mathrm{C}}\approx \mathrm{4}$) by plotting the surface with the operator $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$. The resulting oriented surface precesses about the x axis with a period ${\mathit{\tau }}_{\mathrm{1}}=\mathrm{1}/\stackrel{\mathrm{‾}}{\mathit{\nu }}$. This motion corresponds to the leftmost peak in Fig. 5b.

Now let us consider the action of an antisymmetric part of the measurement operator. The intersection with the z axis of the surface plotted in such a way directly corresponds to the observed J doublet (Fig. 5b). The surface precesses about the x axis with a frequency $\stackrel{\mathrm{‾}}{\mathit{\nu }}$, and its size oscillates with a frequency J (check the corresponding movies to assess the actual evolution of the surface). Such a rotating surface corresponds to a doublet in the ZULF NMR spectrum split around the J frequency (Fig. 5b).

Figure 6 Visualizations of spin dynamics in an AX system (¹H–¹³C nuclear pair) during the ZULF experiment. Conditions used for simulation are the same as that used in Fig. 5 except for assuming an initial density operator of (a) ${\widehat{\mathit{\rho }}}_{\mathrm{0}}=(\mathrm{3}/\mathrm{20})({\widehat{I}}_{\mathrm{2}z}-{\widehat{I}}_{\mathrm{1}z})$ and (b) ${\widehat{\mathit{\rho }}}_{\mathrm{0}}=\frac{\mathrm{1}}{\mathrm{4}}\text{𝟙}+\frac{\mathrm{1}}{\mathrm{4}}({\widehat{I}}_{\mathrm{2}z}+{\widehat{I}}_{\mathrm{1}z})$. Surfaces in panels (a) and (b) are plotted with the measurement operator ${\left({\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$ over a timescale τ₁. The insets of panels (a) and (b) show the evolution of the surfaces over a shorter timescale τ₂. Note the different scaling in panels (a) and (b).

In order to understand the asymmetry of the doublet (Fig. 5b), it is worth analyzing separately the contributions of the symmetric and antisymmetric parts of the initial density matrix to the surface plotted in Fig. 5d. This is now shown in Fig. 6a and b, which visualize the evolution of the antisymmetric, ${\widehat{\mathit{\rho }}}_{\mathrm{0}}=(\mathrm{3}/\mathrm{20})({\widehat{I}}_{\mathrm{2}z}-{\widehat{I}}_{\mathrm{1}z})$, and symmetric, ${\widehat{\mathit{\rho }}}_{\mathrm{0}}=(\mathrm{1}/\mathrm{4})\text{𝟙}+(\mathrm{1}/\mathrm{4})({\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z})$, parts of the initial density operator, respectively. One can see that the evolution of the asymmetric part of the density operator predominantly contributes to the observed ZULF J spectrum (i.e., Fig. 6a is almost equivalent to Fig. 5d); i.e., it gives rise to a symmetric doublet centered at J.

To evaluate this in a more quantitative manner, one can check the intersection of the surface (Fig. 6a) with the z axis. Given the fact that surfaces plotted with the measurement operator ${\left({\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$ are rank-1 spherical harmonics (Appendix D), their shapes can be quantified as p orbitals. Considering the evolution of the surface – i.e., the size being proportional to cos (2πJt) and the angle between the p orbital and the z axis being equal to $(\mathit{\pi }+\mathrm{2}\mathit{\pi }\stackrel{\mathrm{‾}}{\mathit{\nu }}t)$ – the intersection of the surface with the z axis is found to have the following time dependence, $\cos \left(\mathrm{2}\mathit{\pi }Jt\right)\cos (\mathit{\pi }+\mathrm{2}\mathit{\pi }\stackrel{\mathrm{‾}}{\mathit{\nu }}t)=-(\mathrm{1}/\mathrm{2})\left[\cos \left(\mathrm{2}\mathit{\pi }(J+\stackrel{\mathrm{‾}}{\mathit{\nu }})t\right)+\cos \left(\mathrm{2}\mathit{\pi }(J-\stackrel{\mathrm{‾}}{\mathit{\nu }})t\right)\right]$, which indeed gives a symmetric doublet centered at J. Similarly, one can check the contribution of the symmetric part of the density matrix (Fig. 6b) to the measured ZULF J spectrum. In contrast to the surface plotted in Fig. 6a, the surface plotted in Fig. 6b appears $\mathit{\pi }/\mathrm{2}$ radian out of phase in the orientation and $-\mathit{\pi }/\mathrm{2}$ radian out of phase in the size oscillation. This means that the intersection contributing to the measurement has a time dependence $\cos (-\mathit{\pi }/\mathrm{2}+\mathrm{2}\mathit{\pi }Jt)\cos (\mathrm{3}\mathit{\pi }/\mathrm{2}+\mathrm{2}\mathit{\pi }\stackrel{\mathrm{‾}}{\mathit{\nu }}t)=-(\mathrm{1}/\mathrm{2})\left[\cos \left(\mathrm{2}\mathit{\pi }(J+\stackrel{\mathrm{‾}}{\mathit{\nu }})t\right)-\cos \left(\mathrm{2}\mathit{\pi }(J-\stackrel{\mathrm{‾}}{\mathit{\nu }})t\right)\right]$. This term accounts for the small asymmetry of the doublet centered at J shown in Fig. 5b.

Finally, we note that one also needs the surfaces plotted with the measurement operator ${\left({\widehat{\mathrm{ZQ}}}_{\mathit{\pi }/\mathrm{2},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$ (Fig. H1) to completely represent the spin system. Since they do not contribute to the ZULF signal measured with one detector, they are not considered here.

3.3 Signal amplification by reversible exchange

As a second example, we visualize dynamics of the nuclear spin states of dissolved hydrogen during signal amplification by reversible exchange (SABRE) experiments (Adams et al., 2009) at a high magnetic field (Barskiy et al., 2014). The interplay between chemical exchange and coherent spin dynamics is known to induce singlet–triplet mixing (Kiryutin et al., 2017; Knecht et al., 2019; Markelov et al., 2021). For the two protons in free hydrogen (i.e., molecular hydrogen gas dissolved in solution), the singlet–triplet mixing is symmetry-forbidden. However, as soon as the hydrogen molecule binds transiently to a SABRE catalyst and the two protons occupy non-equivalent positions such that chemical symmetry is broken, the difference in their Larmor frequencies, $\mathrm{\Delta }=({\mathit{\delta }}_{\mathrm{1}}-{\mathit{\delta }}_{\mathrm{2}})\mathit{\gamma }{B}_{\mathrm{0}}$, gives rise to singlet–triplet mixing. Such mixing is selective, and only $|S\rangle \to |T_{\mathrm{0}}\rangle$ transitions occur, while the other transitions $|S\rangle \to |T_{\pm }\rangle$ are forbidden at a high field. In the following, we also consider relaxation via locally correlated noise fields to account for the equilibration of triplet sublevels (Ivanov et al., 2008).

Figure 7 Visualizations of spin dynamics in proton pairs of free molecular hydrogen in a high-field SABRE experiment. (a) Scheme of the SABRE experiment: molecular hydrogen initiated in the singlet state undergoes reversible chemical coordination which transiently breaks chemical equivalence between the two spins; NMR, exchange, and relaxation parameters are shown in the inset. (b) Evolution of the spin state populations during the SABRE experiment. (c–f) Visualization of the evolution under the SABRE experiment plotted with the measurement (c) $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$, (d) the singlet population, and (e) ${\left({\widehat{\mathrm{ZQ}}}_{\frac{\mathit{\pi }}{\mathrm{2}},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$. Note the different scaling in panels (c) and (d).

We apply a model recently introduced to simulate spin dynamics in SARBE experiments. All the parameters used for the simulation are given in Fig. 7a, where the chemical rate constants k_a and k_d describe the rates of association and dissociation of H₂, respectively, and the longitudinal relaxation time of the two protons is given by $T_{\mathrm{1}}^{\mathrm{f}}$, when they are free in solution, and $T_{\mathrm{1}}^{\mathrm{b}}$, when they are bound to the SABRE catalyst, respectively. The simulated state populations are shown in Fig. 7b. First, after the start of the exchange process, chemical asymmetry induces the $|S\rangle \to |T_{\mathrm{0}}\rangle$ transition and reduces the population difference between the two states. As the |T₀〉 state is populated, relaxation between triplet sublevels becomes more noticeable, and the states $|T_{\pm }\rangle$ also start to be populated. Since the local noise fields are not perfectly correlated, the populations are finally equalized as shown in Fig. 7b.

The evolving distribution of populations over various states is seen in Fig. 7c–d. In Fig. 7c, the AMPS for the triplet block of the density matrix first grow into an oblate spheroid, which indicates that initially only the |T₀〉 state gets overpopulated. Later, the AMPS evolve into a sphere, indicating that all three triplet states are equally populated as a result of relaxation. Figure 7d shows the changing singlet population. Lastly, Fig. 7e, which measures the out-of-phase coherence ${\left({\widehat{\mathrm{ZQ}}}_{\frac{\mathit{\pi }}{\mathrm{2}},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$, reveals extra information not covered by Fig. 7b; i.e., the singlet–triplet coherence, ${\widehat{I}}_{\mathrm{1}y}{\widehat{I}}_{\mathrm{2}x}-{\widehat{I}}_{\mathrm{1}x}{\widehat{I}}_{\mathrm{2}y}$, is transiently formed during the experiment. Since the J coupling between the two protons is ignored, there is no in-phase coherence (i.e., magnetization differences between the two protons along any direction), and the related surface plot is therefore not included here.

3.4 Spin-lock-induced crossing

Figure 8 shows a typical spin-lock-induced crossing (SLIC) experiment where a continuous radio frequency (RF) field is applied to a pair of chemically inequivalent nuclear spins at a high field, initially prepared in the singlet spin state. The target spins undergo singlet–triplet conversion if the spins are strongly coupled to each other and the amplitude of the RF field matches the J-coupling value between them (DeVience et al., 2013). The parameters used for the simulation and visualization are shown in Fig. 8a, where the molecule is specified by the J-coupling strength between the two spins and Δ, the difference between their Larmor frequencies, expressed in hertz. Under the above conditions, the singlet and triplet populations interconvert with a period of $T=\sqrt{\mathrm{2}}/\mathrm{\Delta }$.

Figure 8 Visualization of spin dynamics in proton pairs during high-field SLIC experiments over a period $T=\sqrt{\mathrm{2}}/\mathrm{\Delta }$. (a) Scheme of the SLIC experiment: an RF field of the amplitude B₁ equal to J coupling between the spins is applied along the rotating-frame x axis at the average proton resonance frequency. (b) Evolution of various spin orders during the SLIC experiment. (c–f) Visualization of the evolution under the SLIC experiment plotted with the measurement (c) $\left(\right|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|{)}_{\widehat{n}}$; (d) the singlet population; (e) ${\left({\widehat{\mathrm{ZQ}}}_{\mathrm{0},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$; (f) ${\left({\widehat{\mathrm{ZQ}}}_{\frac{\mathit{\pi }}{\mathrm{2}},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$. Note the different scaling in panels (c)–(f).

We simulate the SLIC experiment by considering only the coherent processes (no relaxation is included). One could process the simulation results by plotting the population of several chosen operators as a function of time, and the derived population plot is shown in Fig. 8b (note that all operators are defined in the frame rotating in sync with the RF field). Alternatively, one could observe the dynamics and symmetry of the surfaces visualized via our approach to assess the formation of magnetization and/or alignment states. Figure 8c illustrates that the singlet spin order is converted – upon application of a SLIC pulse – into the magnetization opposite to the direction of the RF-field amplitude in the rotating frame. Note here that the conversion of the singlet spin order into magnetization parallel to the direction of the RF field is symmetry-allowed but energy-forbidden, which introduces small but fast oscillations observable through the visualized surfaces (Fig. 8c–f and the movies). Note the formation of in-phase and out-of-phase coherences during application of the SLIC pulse.

3.5 Comparison of the generalized measurement-based visualization and DROPS approach

One should note that the earlier DROPS visualization approach (Garon et al., 2015; Leiner et al., 2017; Leiner and Glaser, 2018) may also be applied to represent the density operator of strongly coupled nuclear spin systems. For DROPS-based visualization, similar to the present approach, one needs to first decompose the density operator into blocks according to the total angular momentum quantum number F. Then each block (F,K) of the density matrix is expanded in terms of the tensor operator basis (see Eq. B1, according to Eq. C1). Next, one maps the expansion of the density operator into a surface expanded over spherical harmonics such that the function of the surface representing the block (F,K) is

$$\begin{array}{}\text{(7)}& {}^{(F,K)}r(\mathit{\theta },\mathit{\varphi })=\sum _{\mathit{\lambda },\mathit{\mu }}{\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}{Y}_{\mathit{\lambda }\mathit{\mu }}(\mathit{\theta },\mathit{\varphi }),\end{array}$$

where ${\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}$ is the polarization moment and Y_λμ(θ,ϕ) is the spherical harmonic. The defined radius can be complex. The amplitudes of ${}^{(F,K)}r(\mathit{\theta },\mathit{\varphi })$ are mapped on the radius into the distance of the surface point along the direction specified by (θ,ϕ), and the phase is mapped into the color of the point according to the color wheel shown in Fig. 9b.

Figure 9 Visualizations of spin dynamics in an AX system (¹H–¹³C nuclear pair) during the ZULF experiment using the multipole tensor operator-based DROPS approach. Conditions are the same as those used in Fig. 5. Panels (a)–(c) show the surfaces representing the spin evolution during the ZULF experiment of (a) the (1,1) block and (b) the (1,0) block. The inset shows the evolution of the surface plotted in panel (b) over a shorter timescale of τ₂. Colors of surfaces denote the phase of the calculated radius (Eq. 7) according to the color wheel shown at the bottom right.

As discussed in the Results section, for a system of N-coupled spins-1/2, the density matrix written in a total angular momentum basis should have in total h² blocks: $h=C_N^{N/\mathrm{2}}$ (binomial coefficient) when N is even and $h=C_N^{(N-\mathrm{1})/\mathrm{2}}$ when N is odd. Given the hermiticity of the density operator, only one of blocks (F,K) or (K,F) needs to be visualized in the DROPS approach. Therefore, only $h(h+\mathrm{1})/\mathrm{2}$ surfaces are needed to fully represent the density operator.

To compare the DROPS approach with the generalized measurement-based visualization approach, we look at the spin dynamics in an AX system during the ZULF experiment once again (Fig. 9). Figure 9a shows the evolution of the surface representing the (1,1) block, which precesses about the x axis with the frequency $\stackrel{\mathrm{‾}}{\mathit{\nu }}$. It resembles the AMPS shown in Fig. 5c apart from a small “blue tail” observed here. Modeling shows this tail to be related to the rank-2 component of the density operator. The (1,0) block is represented through the surface shown in Fig. 9b. The color of the surface changes periodically with the frequency J, and the orientation precesses about the x axis with the frequency $\stackrel{\mathrm{‾}}{\mathit{\nu }}$. While the surface also captures the full dynamics and directly reflects the two characteristic frequencies, J and $\stackrel{\mathrm{‾}}{\mathit{\nu }}$, of the system, it is less obvious how to quantify the measured ZULF J spectra from it.

4 Conclusions

In this paper, we extend an angular momentum probability surfaces (AMPS) approach for visualizing dynamics of quantum spin systems based on generalized observable operators. The plotted 3D surfaces conveniently represent symmetries of density matrices and allow spotting of their presence (orientation, alignment, etc.) or absence even when direct analysis of (time-dependent) density-matrix elements is not obvious. Three different experiments are used to demonstrate the applicability of the novel visualization approach: (i) evolution of the heteronuclear ¹H–¹³C spin pair during the zero- to ultralow-field (ZULF) NMR experiment; (ii) physicochemical conversion of parahydrogen during the signal amplification by reversible exchange (SABRE) experiment at a high magnetic field; (iii) spin dynamics in the ensembles of pairs of spin-1/2 nuclei during the spin-lock-induced crossing (SLIC) experiment. The presented approach allows visualization of complex dynamics in multi-spin systems and may find applications for describing hyperpolarization experiments utilizing parahydrogen (PHIP and SABRE) and general NMR experiments, especially under zero- to ultralow-field conditions.

Appendix A: Proof of rotational invariance of ${\widehat{O}}_{\widehat{n}}$ around $\widehat{n}$

Let us assume that a rotation about the $\widehat{z}$ axis is applied to the measurement observable ${\widehat{O}}_{\widehat{z}}$. Since the rotation keeps the direction $\widehat{z}$ unchanged, according to Eq. (2) the rotated observable can also be used for calculating $r_{\widehat{z}}$. As a result, the two operators for measuring $r_{\widehat{z}}$ must be the same; otherwise, there would be ambiguity for the definition of $r_{\widehat{z}}$.

Similarly, the observable ${\widehat{O}}_{\widehat{n}}$ must be invariant under rotations about $\widehat{n}$ (denoted ${\widehat{R}}_{\widehat{n}}$) to avoid the ambiguity of $r_{\widehat{n}}$. To see this, one can check the following equality:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\widehat{R}}_{\widehat{n}}{\widehat{O}}_{\widehat{n}}{\widehat{R}}_{\widehat{n}}^{-\mathrm{1}}=& {\widehat{R}}_{\widehat{n}}{\widehat{R}}_{\widehat{z}\to \widehat{n}}{\widehat{O}}_{\widehat{z}}{\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}}{\widehat{R}}_{\widehat{n}}^{-\mathrm{1}}\\ =& {\widehat{R}}_{\widehat{z}\to \widehat{n}}\left({\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}}{\widehat{R}}_{\widehat{n}}{\widehat{R}}_{\widehat{z}\to \widehat{n}}\right)\\ & {\widehat{O}}_{\widehat{z}}\left({\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}}{\widehat{R}}_{\widehat{n}}^{-\mathrm{1}}{\widehat{R}}_{\widehat{z}\to \widehat{n}}\right){\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}},\end{array}\end{array}$$

where we have applied Eq. (2) and inserted the operator ${\widehat{R}}_{\widehat{z}\to \widehat{n}}{\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}}$, which equals the identity on both sides.

Given ${\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}}{\widehat{R}}_{\widehat{n}}{\widehat{R}}_{\widehat{z}\to \widehat{n}}={\widehat{R}}_{\widehat{z}}$, one can show that

$$\begin{array}{}\text{(A2)}& {\widehat{R}}_{\widehat{n}}{\widehat{O}}_{\widehat{n}}{\widehat{R}}_{\widehat{n}}^{-\mathrm{1}}={\widehat{R}}_{\widehat{z}\to \widehat{n}}{\widehat{R}}_{\widehat{z}}{\widehat{O}}_{\widehat{z}}{\widehat{R}}_{\widehat{z}}^{-\mathrm{1}}{\widehat{R}}_{\widehat{z}\to \widehat{n}}^{-\mathrm{1}}={\widehat{O}}_{\widehat{n}},\end{array}$$

which proves the conclusion.

Appendix B: Spherical tensor operator basis

B1 Definition of the spherical tensor basis operators

For each block (F,K), a basis of tensor operators is defined by

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}=& \sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\lambda }+\mathrm{1}}{\mathrm{2}F+\mathrm{1}}}}\\ & \sum _{m_F,m_K}\langle K{m}_K\mathit{\lambda }\mathit{\mu }|F{m}_F\rangle |F,m_F\rangle \langle K,m_K|,\end{array}\end{array}$$

with $|F-K|\le \mathit{\lambda }\le F+K$ and $-\mathit{\lambda }\le \mathit{\mu }\le \mathit{\lambda }$ and m_K, m_F, and μ being the projection quantum numbers of the Clebsch–Gordan coefficient 〈Km_Kλμ|Fm_F〉. This derivation generalizes a similar treatment of the diagonal block which can be found in Auzinsh et al. (2010) and Omont (1977).

B2 Properties of the spherical tensor basis operators

The defined operator basis satisfies the following properties.

1. Orthogonality $\mathrm{Tr}\left({\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}{\left({\widehat{\mathbb{T}}}_{{\mathit{\lambda }}^{\prime }{\mathit{\mu }}^{\prime }}^{(F^{\prime },K^{\prime })}\right)}^{†}\right)={\mathit{\delta }}_{FF\prime }{\mathit{\delta }}_{K{K}^{\prime }}{\mathit{\delta }}_{\mathit{\lambda }{\mathit{\lambda }}^{\prime }}{\mathit{\delta }}_{\mathit{\mu }{\mathit{\mu }}^{\prime }}$

2. Completeness: the operator basis is complete (in a sense that it can be used to represent any operator).

3. Sphericalness: given the rotation operator $\widehat{R}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma })=e^{-i\mathit{\alpha }{\widehat{F}}_z}{e}^{-i\mathit{\beta }{\widehat{F}}_y}{e}^{-i\mathit{\gamma }{\widehat{F}}_z}$, a tensor operator transforms under rotation as

   $$\begin{array}{}\text{(B2)}& \begin{array}{rl}& \widehat{R}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma }){\widehat{\mathbb{T}}}_{\mathit{\lambda },\mathit{\mu }}^{(F,K)}{\widehat{R}}^{†}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma })=\\ & \sum _{{\mathit{\mu }}^{\prime }}{\widehat{\mathbb{T}}}_{\mathit{\lambda },{\mathit{\mu }}^{\prime }}^{(F,K)}{D}_{{\mathit{\mu }}^{\prime }\mathit{\mu }}^{\mathit{\lambda }}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma }),\end{array}\end{array}$$

   with Wigner D functions $D_{\mathit{\mu }{\mathit{\mu }}^{\prime }}^{\mathit{\lambda }}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma })=\langle \mathit{\lambda }\mathit{\mu }\left|\widehat{R}\right(\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\left)\right|\mathit{\lambda }{\mathit{\mu }}^{\prime }\rangle$.

4. Hermitian adjoint ${\left({\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{\left(\right(F,K\left)\right)}\right)}^{†}=(-\mathrm{1}{)}^{(F-K+\mathit{\mu })}{\widehat{\mathbb{T}}}_{\mathit{\lambda },-\mathit{\mu }}^{\left(\right(F,K\left)\right)}$

We now prove these properties.

B2.1 Proof of orthogonality

It is obvious that operators defined over different blocks are orthogonal. Next we check the orthogonality between operators within the same block (F,K). For simplicity, the superscripts are omitted.

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& \mathrm{Tr}\left({\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}{\widehat{\mathbb{T}}}_{{\mathit{\lambda }}^{\prime }{\mathit{\mu }}^{\prime }}^{†}\right)={\displaystyle \frac{\sqrt{(\mathrm{2}\mathit{\lambda }+\mathrm{1})(\mathrm{2}{\mathit{\lambda }}^{\prime }+\mathrm{1})}}{\mathrm{2}F+\mathrm{1}}}\\ & \sum _{m_F,m_K}\langle K{m}_K\mathit{\lambda }\mathit{\mu }|F{m}_F\rangle \langle K{m}_K{\mathit{\lambda }}^{\prime }{\mathit{\mu }}^{\prime }|F{m}_F\rangle ,\end{array}\end{array}$$

giving

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}\langle K{m}_K\mathit{\lambda }\mathit{\mu }|F{m}_F\rangle =& (-\mathrm{1}{)}^{K-m_K}\sqrt{{\displaystyle \frac{\mathrm{2}F+\mathrm{1}}{\mathrm{2}\mathit{\lambda }+\mathrm{1}}}}\\ & \langle F{m}_{F}K(-m_K)|\mathit{\lambda }\mathit{\mu }\rangle .\end{array}\end{array}$$

The sum can be evaluated as

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}& \sum _{m_K,m_F}\langle K{m}_K\mathit{\lambda }\mathit{\mu }|F{m}_F\rangle \langle K{m}_K{\mathit{\lambda }}^{\prime }{\mathit{\mu }}^{\prime }|F{m}_F\rangle \\ & ={\displaystyle \frac{\mathrm{2}F+\mathrm{1}}{\mathrm{2}\mathit{\lambda }+\mathrm{1}}}\sum _{m_K,m_F}\langle \mathit{\lambda }\mathit{\mu }|F{m}_{F}K{m}_K\rangle \langle F{m}_{F}K{m}_K|{\mathit{\lambda }}^{\prime }{\mathit{\mu }}^{\prime }\rangle \\ & ={\displaystyle \frac{\mathrm{2}F+\mathrm{1}}{\mathrm{2}\mathit{\lambda }+\mathrm{1}}}{\mathit{\delta }}_{\mathit{\lambda }{\mathit{\lambda }}^{\prime }}{\mathit{\delta }}_{\mathit{\mu }{\mathit{\mu }}^{\prime }}.\end{array}\end{array}$$

B2.2 Proof of completeness

Since the operators are orthogonal to each other, they must be linearly independent. The total number of operators over block (F,K) (assuming F>K without lose of generality) is

$$\begin{array}{}\text{(B6)}& \sum _{\mathit{\lambda }=F-K}^{F+K}\mathrm{2}\mathit{\lambda }+\mathrm{1}=(\mathrm{2}F+\mathrm{1})(\mathrm{2}K+\mathrm{1}).\end{array}$$

Since the number of independent basis operators equals the number of degrees of freedom over the block (see Appendix B1), the defined operator basis is complete.

B2.3 Proof of sphericalness

It is equivalent to prove the following conditions (the superscript (F,K) is omitted):

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}& [{\widehat{F}}_z,{\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}]=\mathit{\mu }{\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }},\\ & [{\widehat{F}}_{\pm },{\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}]=\sqrt{\mathit{\lambda }(\mathit{\lambda }+\mathrm{1})-\mathit{\mu }(\mathit{\mu }\pm \mathrm{1})}{\widehat{\mathbb{T}}}_{\mathit{\lambda },\mathit{\mu }\pm \mathrm{1}}.\end{array}\end{array}$$

Proof of the first equation is obvious. For the second part, we have

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}& [{\widehat{F}}_{\pm },{\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}]\\ & =\sum _{m_F,m_K}\langle K{m}_K\mathit{\lambda }\mathit{\mu }|F{m}_F\rangle \\ & \left({\widehat{F}}_{\pm }|F,m_F\rangle \langle K,m_K|-|F,m_F\rangle \langle K,m_K|{\widehat{F}}_{\pm }\right)\\ & =(\sqrt{F(F+\mathrm{1})-m_F(m_F\mp \mathrm{1})}\langle K{m}_K\mathit{\lambda }\mathit{\mu }|F(m_F\mp \mathrm{1})\rangle \\ & -\sqrt{K(K+\mathrm{1})-m_K(m_K\pm \mathrm{1})}\langle K(m_K\pm \mathrm{1})\mathit{\lambda }\mathit{\mu }|F{m}_F\rangle )\\ & |F,m_F\rangle \langle K,m_K|.\end{array}\end{array}$$

We can apply Eq. (B4) and the recursion relationships between Clebsch–Gordan coefficients to find that

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}& C_{\pm }(F,m_F\mp \mathrm{1})\langle F{m}_{F}K(-m_K\mp \mathrm{1})|\mathit{\lambda }\mathit{\mu }\rangle \\ & +C_{\pm }(K,-m_K\mp \mathrm{1})\langle F(m_F\mp \mathrm{1})K(-m_K)|\mathit{\lambda }\mathit{\mu }\rangle \\ & =C_{\pm }(\mathit{\lambda },\mathit{\mu })\langle F{m}_{F}K(-m_K)|\mathit{\lambda }\mathit{\mu }\pm \mathrm{1}\rangle .\end{array}\end{array}$$

With $C_{\pm }(F,m)=\sqrt{F(F+\mathrm{1})-m(m\pm \mathrm{1})}$, Eq. (B8) can be simplified to yield the second part of Eq. (B7).

B2.4 Proof of a Hermitian adjoint

Lastly, we check the Hermitian adjoint property of the basis operators:

$$\begin{array}{}\text{(B10)}& \begin{array}{rl}& {\left({\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\right)}^{†}=\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\lambda }+\mathrm{1}}{\mathrm{2}F+\mathrm{1}}}}\\ & \sum _{m_F,m_K}\langle K{m}_K\mathit{\lambda }\mathit{\mu }|F{m}_F\rangle |K,m_K\rangle \langle F,m_F|,\\ & {\widehat{\mathbb{T}}}_{\mathit{\lambda },-\mathit{\mu }}^{(K,F)}=\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\lambda }+\mathrm{1}}{\mathrm{2}K+\mathrm{1}}}}\\ & \sum _{m_F,m_K}\langle F{m}_F\mathit{\lambda }(-\mathit{\mu })|K{m}_K\rangle |K,m_K\rangle \langle F,m_F|.\end{array}\end{array}$$

Given Eq. (B4) and the symmetry of the Clebsch–Gordan coefficients

$$\begin{array}{}\text{(B11)}& \langle F{m}_{F}K(-m_K)|\mathit{\lambda }\mathit{\mu }\rangle =\langle K{m}_{K}F(-m_F\left)\right|\mathit{\lambda }(-\mathit{\mu })\rangle ,\end{array}$$

we have

$$\begin{array}{}\text{(B12)}& {\left({\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\right)}^{†}=(-\mathrm{1}{)}^{(F-K+\mathit{\mu })}{\widehat{\mathbb{T}}}_{\mathit{\lambda },-\mathit{\mu }}^{(K,F)}.\end{array}$$

Appendix C: Expansion of the density operator over the spherical tensor operator basis

The density operator could be expanded over the spherical tensor operator basis:

$$\begin{array}{}\text{(C1)}& \widehat{\mathit{\rho }}=\sum _{F,K}\sum _{\mathit{\lambda }=|F-K|}^{F+K}\sum _{\mathit{\mu }=-\mathit{\lambda }}^{\mathit{\lambda }}{\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}{\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}.\end{array}$$

Applying the orthogonality condition of the basis operators (Appendix B2), we have the so-called polarization moments (Auzinsh et al., 2010):

$$\begin{array}{}\text{(C2)}& {\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}=\mathrm{Tr}\left[\widehat{\mathit{\rho }}{\left({\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\right)}^{†}\right].\end{array}$$

Appendix D: Completeness of the visualization

The visualization is complete in a sense that there is a one-to-one correspondence between the surfaces and the related density operators. First, as discussed in the main text, the surfaces can be uniquely determined from the given density operator. Here we show that the visualized density operator can be reconstructed from the surfaces. Following Eqs. (5) and (6) and formulating the rotating operator ${\widehat{R}}_{\widehat{z}\to \widehat{n}}$ as $\widehat{R}(\mathit{\varphi },\mathit{\theta },\mathrm{0})$ (see Appendix B2 for a definition and assume $\widehat{n}$ is pointing at (θ,ϕ)), the surfaces can be expressed using the function $f_m^{(F,K)}(\mathit{\theta },\mathit{\varphi })=\mathrm{Tr}\left(\widehat{\mathit{\rho }}\widehat{R}(\mathit{\varphi },\mathit{\theta },\mathrm{0})|K,m\rangle \langle F,m|{\widehat{R}}^{-\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})\right)$ as

$$\begin{array}{}\text{(D1)}& \begin{array}{rl}& {}_{\mathit{\phi },m}^{(F,K)}r(\mathit{\theta },\mathit{\varphi })\\ & ={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\left(e^{i\mathit{\phi }}{f}_m^{(F,K)}(\mathit{\theta },\mathit{\varphi })+e^{-i\mathit{\phi }}{f}_m^{(K,F)}(\mathit{\theta },\mathit{\varphi })\right).\end{array}\end{array}$$

There are two useful proprieties of the defined function $f_m^{(F,K)}$.

1. Expansion over spherical harmonics

   $$\begin{array}{}\text{(D2)}& \begin{array}{rl}{f}_m^{(F,K)}=& (-\mathrm{1}{)}^{Fk-m}\sum _{\mathit{\lambda },\mathit{\mu }}\langle FmK(-m\left)\right|\mathit{\lambda }\mathrm{0}\rangle \\ & \sqrt{{\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathrm{2}\mathit{\lambda }+\mathrm{1}}}}{\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}{Y}_{\mathit{\lambda }\mathit{\mu }},\end{array}\end{array}$$

   where ${\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}$ are the polarization moments (Auzinsh et al., 2010) and Y_λμ(θ,ϕ) are spherical harmonics (Garon et al., 2015).

2. Complex conjugation

   $$\begin{array}{}\text{(D3)}& f_m^{(K,F)}(\mathit{\theta },\mathit{\varphi })={\left(f_m^{(F,K)}(\mathit{\theta },\mathit{\varphi })\right)}^{*}\end{array}$$

The first propriety (see the proof below) allows us to extract the polarization moments of the density operator from the function $f_m^{(F,K)}$ by decomposing it into spherical harmonics and hence the density operator itself (see Appendix C). Function $f_m^{(F,K)}$ is calculated from the function of the surfaces plotted.

For a diagonal block (F,F), the function $f_m^{(F,F)}$ is the same as the function of AMPS by setting $\mathit{\phi }=\mathit{\pi }/\mathrm{4}$ and m=F:

$$\begin{array}{}\text{(D4)}& f_F^{(F,F)}{=}_{\mathit{\pi }/\mathrm{4},F}^{(F,F)}r.\end{array}$$

For an off-diagonal block (F,K), when two surfaces are plotted with two measurement observables ${\widehat{O}}_{\mathrm{0},m}^{(F,K)}$ and ${\widehat{O}}_{\mathit{\pi }/\mathrm{2},m}^{(F,K)}$, one can derive the function $f_m^{(F,K)}$ by

$$\begin{array}{}\text{(D5)}& f_m^{(F,K)}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\left({}_{\mathrm{0},m}^{(F,K)}r-i{}_{\mathit{\pi }/\mathrm{2},m}^{(F,K)}r\right).\end{array}$$

Also, we note here that the second propriety guarantees that the surfaces plotted will always be real.

D1 Proof of the expansion of the function $f_m^{(F,K)}$ over spherical harmonics

Substituting the expansion of the density operator Eq. (C1) into the definition of the function $f_m^{(F,K)}(\mathit{\theta },\mathit{\varphi })$ gives

$$\begin{array}{}\text{(D6)}& \begin{array}{rl}& f_m^{(F,K)}(\mathit{\theta },\mathit{\varphi })\\ & =\sum _{\mathit{\lambda },\mathit{\mu }}{\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\langle F,m\left|{\widehat{R}}^{-\mathrm{1}}\right(\mathit{\varphi },\mathit{\theta },\mathrm{0}\left){\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\widehat{R}\right(\mathit{\varphi },\mathit{\theta },\mathrm{0}\left)\right|K,m\rangle \\ & =\sum _{\mathit{\lambda },\mathit{\mu }}{\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\sum _{{\mathit{\mu }}^{\prime }}\langle F,m\left|{\widehat{\mathbb{T}}}_{\mathit{\lambda }{\mathit{\mu }}^{\prime }}^{(F,K)}\right|K,m\rangle D_{{\mathit{\mu }}^{\prime }\mathit{\mu }}^{\mathit{\lambda }}(\mathrm{0},-\mathit{\theta },-\mathit{\varphi }),\end{array}\end{array}$$

where we have employed the sphericalness of the basis operators (Appendix B2). The term $\langle F,m\left|{\widehat{\mathbb{T}}}_{\mathit{\lambda }{\mathit{\mu }}^{\prime }}^{(F,K)}\right|K,m\rangle$ can be easily estimated from the definition (B1). Note that μ^′ can only take a zero value. In this particular case, Wigner D functions are related to spherical harmonics via

$$\begin{array}{}\text{(D7)}& \begin{array}{rl}{D}_{\mathrm{0}\mathit{\mu }}^{\mathit{\lambda }}(\mathrm{0},-\mathit{\theta },-\mathit{\varphi })& ={\left(D_{\mathit{\mu }\mathrm{0}}^{\mathit{\lambda }}(\mathit{\varphi },\mathit{\theta },\mathrm{0})\right)}^{*}\\ & =\sqrt{{\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathrm{2}\mathit{\lambda }+\mathrm{1}}}}{Y}_{\mathit{\lambda }\mathit{\mu }}(\mathit{\theta },\mathit{\varphi }).\end{array}\end{array}$$

Substituting Eq. (D7) into Eq. (D6) proves Eq. (D2).

D2 Proof of the complex conjugation propriety of the function $f_m^{(F,K)}$

Starting from Eq. (D2), by exchanging K and F, one can find that

$$\begin{array}{}\text{(D8)}& \begin{array}{rl}{f}_m^{(K,F)}=& (-\mathrm{1}{)}^{Fl-m}\sum _{\mathit{\lambda },\mathit{\mu }}\langle KmF(-m\left)\right|\mathit{\lambda }\mathrm{0}\rangle \\ & \sqrt{{\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathrm{2}\mathit{\lambda }+\mathrm{1}}}}{\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(K,F)}{Y}_{\mathit{\lambda }\mathit{\mu }}(\mathit{\theta },\mathit{\varphi }).\end{array}\end{array}$$

To find out the relationship between ${\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(K,F)}$ and ${\mathit{\rho }}_{\mathit{\lambda },\mathit{\mu }}^{(F,K)}$, we employ Eq. (B12):

$$\begin{array}{}\text{(D9)}& \begin{array}{rl}{\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(K,F)}& =\mathrm{Tr}\left(\widehat{\mathit{\rho }}{\left({\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(K,F)}\right)}^{†}\right)\\ & =(-\mathrm{1}{)}^{(K-F+\mathit{\mu })}\mathrm{Tr}\left(\widehat{\mathit{\rho }}{\widehat{\mathbb{T}}}_{\mathit{\lambda },-\mathit{\mu }}^{(F,K)}\right)\\ & =(-\mathrm{1}{)}^{(K-F+\mathit{\mu })}{\left(\mathrm{Tr}\left(\widehat{\mathit{\rho }}{\left({\widehat{\mathbb{T}}}_{\mathit{\lambda },-\mathit{\mu }}^{(F,K)}\right)}^{†}\right)\right)}^{*}\\ & =(-\mathrm{1}{)}^{(K-F+\mathit{\mu })}{\left({\mathit{\rho }}_{\mathit{\lambda },-\mathit{\mu }}^{(F,K)}\right)}^{*}.\end{array}\end{array}$$

By inserting

$$\begin{array}{}\text{(D10)}& {\left(Y_{\mathit{\lambda }\mathit{\mu }}(\mathit{\theta },\mathit{\varphi })\right)}^{*}=(-\mathrm{1}{)}^{\mathit{\mu }}{Y}_{\mathit{\lambda },-\mathit{\mu }}(\mathit{\theta },\mathit{\varphi })\end{array}$$

and

$$\begin{array}{}\text{(D11)}& \langle KmF(-m)|\mathit{\lambda }\mathrm{0}\rangle =\langle FmK(-m\left)\right|\mathit{\lambda }\mathrm{0}\rangle \end{array}$$

into Eq. (D8) and substituting the summation label μ with −μ, one can derive Eq. (D3).

Appendix E: Rotational symmetry

Certain rotational symmetries of the density operator can be directly reflected from the visualized surfaces.

1. Global rotation. Any global rotation can be directly reflected by the rotation of the surfaces.

2. Higher-order symmetries. All the visualized surfaces have a q-fold symmetry around $\widehat{z}$, and then all the elements of ${\widehat{\mathit{\rho }}}_{m_F,m_K}^{(F,K)}$ are to be zero, except those for which $m_F-m_K=qN$, where N is an integer.

E1 Proof of the propriety of surfaces under global rotations

Assume a global rotation is applied; the resulting rotated density operator is formulated as

$$\begin{array}{}\text{(E1)}& {\widehat{\mathit{\rho }}}^{\prime }=\widehat{R}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma })\widehat{\mathit{\rho }}{\widehat{R}}^{†}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma }).\end{array}$$

Following from the sphericalness (Appendix B2) of the basis operators, the polarization moments of ${\widehat{\mathit{\rho }}}^{\prime }$ are

$$\begin{array}{}\text{(E2)}& {\left({\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\right)}^{\prime }=D_{\mathit{\mu }{\mathit{\mu }}^{\prime }}^{\mathit{\lambda }}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma }){\mathit{\rho }}_{\mathit{\lambda }{\mathit{\mu }}^{\prime }}^{(F,K)}.\end{array}$$

Inserting ${\left({\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\right)}^{\prime }$ into Eq. (D2), one can derive the function ${\left(f_m^{(F,K)}\right)}^{\prime }$ defined with respect to the rotated density operator ${\widehat{\mathit{\rho }}}^{\prime }$:

$$\begin{array}{}\text{(E3)}& \begin{array}{rl}{\left(f_m^{(F,K)}\right)}^{\prime }=& (-\mathrm{1}{)}^{Fk-m}\sum _{\mathit{\lambda },\mathit{\mu }}\langle FmK(-m\left)\right|\mathit{\lambda }\mathrm{0}\rangle \\ & \sqrt{{\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathrm{2}\mathit{\lambda }+\mathrm{1}}}}{\left({\mathit{\rho }}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\right)}^{\prime }{Y}_{\mathit{\lambda }\mathit{\mu }}(\mathit{\theta },\mathit{\varphi }).\end{array}\end{array}$$

On the other hand, one can apply the same rotation to surface $\left(f_m^{(F,K)}\right)$ and derive a new surface ${\left(f_m^{(F,K)}\right)}^{\prime \prime }$:

$$\begin{array}{}\text{(E4)}& {\left(f_m^{(F,K)}\right)}^{\prime \prime }=\widehat{R}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma })f_m^{(F,K)}.\end{array}$$

Following from the rotational symmetry of spherical harmonics,

$$\begin{array}{}\text{(E5)}& \begin{array}{rl}\widehat{R}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma })Y_{\mathit{\lambda }\mathit{\mu }}(\mathit{\theta },\mathit{\varphi })& =Y_{\mathit{\lambda }\mathit{\mu }}\left({\widehat{R}}^{-\mathrm{1}}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma })(\mathit{\theta },\mathit{\varphi })\right)\\ & =\sum _{{\mathit{\mu }}^{\prime }}{Y}_{\mathit{\lambda }{\mathit{\mu }}^{\prime }}(\mathit{\theta },\mathit{\varphi })D_{{\mathit{\mu }}^{\prime }\mathit{\mu }}^{\mathit{\lambda }}(\mathit{\alpha },\mathit{\beta },\mathit{\gamma }),\end{array}\end{array}$$

one can realize ${\left(f_m^{(F,K)}\right)}^{\prime \prime }={\left(f_m^{(F,K)}\right)}^{\prime }$, which proves the conclusion.

E2 Proof of the connection between density-matrix coherence and higher-order symmetries of the surfaces

If a q-fold symmetry is presented for all the surfaces, the density operator would also have a q-fold symmetry since they rotate in sync, i.e.,

$$\begin{array}{}\text{(E6)}& {\widehat{\mathit{\rho }}}_{m_F,m_K}^{(F,K)}{e}^{-i(m_F-m_K)\frac{\mathrm{2}\mathit{\pi }}{q}}={\widehat{\mathit{\rho }}}_{m_F,m_K}^{(F,K)}.\end{array}$$

This requires all the elements of ${\widehat{\mathit{\rho }}}_{m_F,m_K}^{(F,K)}$ to be zero except those for which $m_F-m_K=qN$, where N is an integer.

Appendix F: Step-by-step implementation

Here we give step-by-step instructions on plotting the AMPS for the density operator $\widehat{\mathit{\rho }}=\frac{\mathrm{1}}{\mathrm{4}}\text{𝟙}+\frac{\mathrm{1}}{\mathrm{4}}({\widehat{I}}_{\mathrm{1}z}+{\widehat{I}}_{\mathrm{2}z})$.

Step 1: find out the measurement operator along the $\widehat{z}$ axis ${\widehat{O}}_{\widehat{z}}$. For the example shown here, ${\widehat{O}}_{\widehat{z}}=|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|$.

Step 2: write the operator ${\widehat{O}}_{\widehat{z}}$ in the total angular momentum basis ($|\mathrm{1},\mathrm{1}\rangle$, $|\mathrm{1},\mathrm{0}\rangle$, $|\mathrm{1},-\mathrm{1}\rangle$, and $|\mathrm{0},\mathrm{0}\rangle$).

$$\begin{array}{}\text{(F1)}& {\widehat{O}}_{\widehat{z}}=|\mathrm{1},\mathrm{1}\rangle \langle \mathrm{1},\mathrm{1}|=\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right).\end{array}$$

Step 3: construct the rotating operator $\widehat{R}(\mathit{\varphi },\mathit{\theta },\mathrm{0})$ (see Appendix B2) in the total angular momentum basis. Only the diagonal blocks (F,F) have nonzero entries, i.e., $\langle F,m_{\mathrm{1}}\left|\widehat{R}\right(\mathit{\varphi },\mathit{\theta },\mathrm{0}\left)\right|F,m_{\mathrm{2}}\rangle =D_{m_{\mathrm{1}},m_{\mathrm{2}}}^F(\mathit{\varphi },\mathit{\theta },\mathrm{0})$. For reference, the Wigner D matrix $D_{m_{\mathrm{1}},m_{\mathrm{2}}}^F(\mathit{\varphi },\mathit{\theta },\mathrm{0})$ can be calculated in Mathematica using the function Wigner D$\left[\mathit{\right\{}F,m_{\mathrm{1}},m_{\mathrm{2}}\mathit{\}},-\mathit{\varphi },-\mathit{\theta },\mathrm{0}]$. For the example here,

$$\begin{array}{}\text{(F2)}& \begin{array}{rl}& \widehat{R}(\mathit{\varphi },\mathit{\theta },\mathrm{0})\\ & =\left(\begin{array}{cccc}{D}_{\mathrm{1},\mathrm{1}}^{\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})& D_{\mathrm{1},\mathrm{0}}^{\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})& D_{\mathrm{1},-\mathrm{1}}^{\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})& \mathrm{0}\\ D_{\mathrm{0},\mathrm{1}}^{\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})& D_{\mathrm{0},\mathrm{0}}^{\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})& D_{\mathrm{0},-\mathrm{1}}^{\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})& \mathrm{0}\\ D_{-\mathrm{1},\mathrm{1}}^{\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})& D_{-\mathrm{1},\mathrm{0}}^{\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})& D_{-\mathrm{1},-\mathrm{1}}^{\mathrm{1}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& D_{\mathrm{0},\mathrm{0}}^{\mathrm{0}}(\mathit{\varphi },\mathit{\theta },\mathrm{0})\end{array}\right)\\ & =\left(\begin{array}{cccc}\frac{\mathrm{1}}{\mathrm{2}}{e}^{-i\mathit{\varphi }}\left(\cos \right(\mathit{\theta })+\mathrm{1})& -\frac{e^{-i\mathit{\varphi }}\sin \left(\mathit{\theta }\right)}{\sqrt{\mathrm{2}}}& e^{-i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)& \mathrm{0}\\ \frac{\sin \left(\mathit{\theta }\right)}{\sqrt{\mathrm{2}}}& \cos \left(\mathit{\theta }\right)& -\frac{\sin \left(\mathit{\theta }\right)}{\sqrt{\mathrm{2}}}& \mathrm{0}\\ e^{i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)& \frac{e^{i\mathit{\varphi }}\sin \left(\mathit{\theta }\right)}{\sqrt{\mathrm{2}}}& \frac{\mathrm{1}}{\mathrm{2}}{e}^{i\mathit{\varphi }}\left(\cos \right(\mathit{\theta })+\mathrm{1})& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right).\end{array}\end{array}$$

Step 4: construct the measurement operator ${\widehat{O}}_{\widehat{n}}$ along $\widehat{n}$ ($\widehat{n}$ points to (θ,ϕ)) using Eq. (2) and ${\widehat{R}}_{\widehat{z}\to \widehat{n}}=\widehat{R}(\mathit{\varphi },\mathit{\theta },\mathrm{0})$ (see Appendix B2 for the definition of $\widehat{R}(\mathit{\varphi },\mathit{\theta },\mathrm{0})$):

$$\begin{array}{}\text{(F3)}& {\widehat{O}}_{\widehat{n}}={\widehat{R}}_{(}\mathit{\varphi },\mathit{\theta },\mathrm{0}\left){\widehat{O}}_{\widehat{z}}{\widehat{R}}^{-\mathrm{1}}\right(\mathit{\varphi },\mathit{\theta },\mathrm{0}).\end{array}$$

The operator ${\widehat{O}}_{\widehat{n}}$ in this example can be calculated using the matrix derived in Eq. (F1) and Eq. (F2) as

$$\begin{array}{}\text{(F4)}& \begin{array}{rl}& {\widehat{O}}_{\widehat{n}}\\ & =\left(\begin{array}{cccc}{\cos }^{\mathrm{4}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)& \frac{e^{-i\mathit{\varphi }}\sin \left(\mathit{\theta }\right)\left(\cos \right(\mathit{\theta })+\mathrm{1})}{\mathrm{2}\sqrt{\mathrm{2}}}& \frac{\mathrm{1}}{\mathrm{4}}{e}^{-\mathrm{2}i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\mathit{\theta }\right)& \mathrm{0}\\ \frac{e^{i\mathit{\varphi }}\sin \left(\mathit{\theta }\right)\left(\cos \right(\mathit{\theta })+\mathrm{1})}{\mathrm{2}\sqrt{\mathrm{2}}}& \frac{{\sin }^{\mathrm{2}}\left(\mathit{\theta }\right)}{\mathrm{2}}& \frac{e^{-i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)\sin \left(\mathit{\theta }\right)}{\sqrt{\mathrm{2}}}& \mathrm{0}\\ \frac{\mathrm{1}}{\mathrm{4}}{e}^{\mathrm{2}i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\mathit{\theta }\right)& \frac{e^{i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)\sin \left(\mathit{\theta }\right)}{\sqrt{\mathrm{2}}}& {\sin }^{\mathrm{4}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right).\end{array}\end{array}$$

Step 5: write the density operator in Zeeman basis (|αα〉, |αβ〉, |βα〉, and |ββ〉).

$$\begin{array}{}\text{(F5)}& {\widehat{\mathit{\rho }}}_{\left(\mathrm{Zeeman}\right)}=\left(\begin{array}{cccc}\frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{\mathrm{4}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{4}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\end{array}$$

Step 6: write the transformation matrix P accounting for the basis switching from Zeeman basis to total angular momentum basis. Generally, the entries P_ij equal 〈ϕ_i|ψ_j〉 when switching the basis vectors, $\mathit{\left\{}\right|{\mathit{\varphi }}_{\mathrm{1}}\rangle ,|{\mathit{\varphi }}_{\mathrm{2}}\rangle ,\mathrm{\dots }\mathit{\}}$, to another set of basis vectors, $\mathit{\left\{}\right|{\mathit{\psi }}_{\mathrm{1}}\rangle ,|{\mathit{\psi }}_{\mathrm{2}}\rangle ,\mathrm{\dots }\mathit{\}}$.

$$\begin{array}{}\text{(F6)}& P=\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}& \mathrm{0}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}\\ \mathrm{0}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}& \mathrm{0}& -\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right)\end{array}$$

Step 7: write the density operator in total angular momentum basis using ${\widehat{\mathit{\rho }}}_{\left(\mathrm{SingletTriplet}\right)}=P^{-\mathrm{1}}{\widehat{\mathit{\rho }}}_{\left(\mathrm{Zeeman}\right)}P$.

$$\begin{array}{}\text{(F7)}& \begin{array}{rl}& {\widehat{\mathit{\rho }}}_{\left(\mathrm{SingletTriplet}\right)}\\ & =\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}& -\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}& \mathrm{0}\end{array}\right)\left(\begin{array}{cccc}\frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{\mathrm{4}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{4}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\\ & \left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}& \mathrm{0}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}\\ \mathrm{0}& \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}& \mathrm{0}& -\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right)\\ & =\left(\begin{array}{cccc}\frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{\mathrm{4}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{4}}\end{array}\right)\end{array}\end{array}$$

Step 8: calculate the distance of the surface to the origin along $\widehat{n}$ using $r(\mathit{\theta },\mathit{\varphi })=\mathrm{Tr}\left({\widehat{O}}_{\widehat{n}}{\widehat{\mathit{\rho }}}_{\left(\mathrm{SingletTriplet}\right)}\right)$. In the simple example given here, the results can be calculated symbolically. Note that, for more complicated cases, this can be done numerically.

$$\begin{array}{}\text{(F8)}& \begin{array}{rl}& r(\mathit{\theta },\mathit{\varphi })\\ & =\mathrm{Tr}\left(\left(\begin{array}{cccc}\frac{\mathrm{1}}{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{\mathrm{4}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{\mathrm{4}}\end{array}\right)\right.\\ & \left.\left(\begin{array}{cccc}{\cos }^{\mathrm{4}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)& \frac{e^{-i\mathit{\varphi }}\sin \left(\mathit{\theta }\right)\left(\cos \right(\mathit{\theta })+\mathrm{1})}{\mathrm{2}\sqrt{\mathrm{2}}}& \frac{\mathrm{1}}{\mathrm{4}}{e}^{-\mathrm{2}i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\mathit{\theta }\right)& \mathrm{0}\\ \frac{e^{i\mathit{\varphi }}\sin \left(\mathit{\theta }\right)\left(\cos \right(\mathit{\theta })+\mathrm{1})}{\mathrm{2}\sqrt{\mathrm{2}}}& \frac{{\sin }^{\mathrm{2}}\left(\mathit{\theta }\right)}{\mathrm{2}}& \frac{e^{-i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)\sin \left(\mathit{\theta }\right)}{\sqrt{\mathrm{2}}}& \mathrm{0}\\ \frac{\mathrm{1}}{\mathrm{4}}{e}^{\mathrm{2}i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\mathit{\theta }\right)& \frac{e^{i\mathit{\varphi }}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)\sin \left(\mathit{\theta }\right)}{\sqrt{\mathrm{2}}}& {\sin }^{\mathrm{4}}\left(\frac{\mathit{\theta }}{\mathrm{2}}\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}(\mathrm{1}+\cos \mathit{\theta })\end{array}\end{array}$$

Step 9: plot the surface using the calculated radius along various directions.

Appendix G: Equivalence between the surface function of AMPS and the Husimi Q function

One may see the equivalence by comparing the measurement observables directly. For a diagonal block (F,F), the measurement observable ${\widehat{O}}_{\widehat{z}}$ for AMPS is $|F,F\rangle \langle F,F|$, which can be expanded over the spherical tensor operator basis (see Appendix B1) as

$$\begin{array}{}\text{(G1)}& |F,F\rangle \langle F,F|=\sum _{\mathit{\lambda },\mathit{\mu }}{c}_{\mathit{\lambda },\mathit{\mu }}{\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(F,F)},\end{array}$$

where the coefficient c_λ,μ can be calculated using Eq. (C2):

$$\begin{array}{}\text{(G2)}& \begin{array}{rl}{c}_{\mathit{\lambda },\mathit{\mu }}& =\mathrm{Tr}\left[|F,F\rangle \langle F,F|{\left({\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathit{\mu }}^{(F,K)}\right)}^{†}\right]\\ & =\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\lambda }+\mathrm{1}}{\mathrm{2}F+\mathrm{1}}}}\langle FF\mathit{\lambda }\mathit{\mu }|FF\rangle .\end{array}\end{array}$$

The Clebsch–Gordan coefficient 〈FFλμ|FF〉 would only be nonzero when $\mathrm{0}\le \mathit{\lambda }\le \mathrm{2}F$. When λ is within this range, one may find

$$\begin{array}{}\text{(G3)}& \begin{array}{rl}& \langle FF\mathit{\lambda }\mathit{\mu }|FF\rangle =\\ & {\mathit{\delta }}_{\mathit{\mu },\mathrm{0}}\sqrt{\mathrm{2}F+\mathrm{1}}\left(\mathrm{2}F\right)\mathrm{!}\left[\right(\mathrm{2}F-\mathit{\lambda }\left)\mathrm{!}\right(\mathrm{1}+\mathrm{2}F+\mathit{\lambda })\mathrm{!}{]}^{-\mathrm{1}/\mathrm{2}}.\end{array}\end{array}$$

Inserting it into Eq. (G2),

$$\begin{array}{}\text{(G4)}& \begin{array}{rl}& |F,F\rangle \langle F,F|=\\ & \sum _{\mathit{\lambda }=\mathrm{0}}^{\mathrm{2}F}\sqrt{\mathrm{2}\mathit{\lambda }+\mathrm{1}}\left(\mathrm{2}F\right)\mathrm{!}\left[\right(\mathrm{2}F-\mathit{\lambda }\left)\mathrm{!}\right(\mathrm{1}+\mathrm{2}F+\mathit{\lambda })\mathrm{!}{]}^{-\mathrm{1}/\mathrm{2}}\\ & {\widehat{\mathbb{T}}}_{\mathit{\lambda }\mathrm{0}}^{(F,F)},\end{array}\end{array}$$

which is the same as the operator used for calculating the Husimi Q function (Koczor et al., 2020, Eq. 5).

Appendix H: Surfaces plotted with ${\left({\widehat{\mathrm{ZQ}}}_{\mathit{\pi }/\mathrm{2},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$ describing the spin dynamics in an AX system during ZULF measurement

Figure H1 Visualizations of spin dynamics in an AX system (¹H–¹³C nuclear pair) during the zero- to ultralow-field (ZULF) nuclear magnetic resonance experiment given the same conditions used in Fig. 5. Surfaces are plotted with the measurement operator ${\left({\widehat{\mathrm{ZQ}}}_{\frac{\mathit{\pi }}{\mathrm{2}},\mathrm{0}}^{(\mathrm{1},\mathrm{0})}\right)}_{\widehat{n}}$ over a timescale τ₁; the inset shows the evolution of the surfaces over a shorter timescale τ₂.