Visualization of dynamics in coupled multi-spin systems

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


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 . 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 = −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 (φ, θ, 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. The visualized surface is a collection of points plotted at a distance rn from the origin along the directionn (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.

Results
In our visualization approach, a density matrixρ is represented by the action of a measurement operatorÔn: rn = Tr(ρÔn). (1) Here, rn indicates the distance of the plotted surface point to the origin alongn (Fig. 1a). In order to find the exact form of the operatorÔn, a measurement is first defined aŝ Oẑ (note that the subscriptẑ denotes the direction of plotting the result of the measurement and does not specify the direction of the measured property). For example,Ôẑ could be (Î 1xÎ2y −Î 1yÎ2x ); see below. The operatorÔn then represents the observable obtained via global rotationRẑ →n , which bringsẑ ton: There is one important caveat. The observableÔẑ must remain unchanged under the rotation aboutẑ; 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Ôn must be Hermitian for a measurement rn to be real. Under such constraints, one may still obtain a negative value for rn; thus, color is introduced as an additional "degree of freedom". As a convention, in this work we set the distance of the surface alongn as the amplitude of rn and set the color of the surface alongn to be red if rn is positive and blue if rn 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 {Ôn} 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 operatorsF z and (F,F) =F 2 x +F 2 y +F 2 z =F 2 : whereF α is a total spin projection operator (α = x, y, z), F α = N i=1Î iα for N spins-1/2, andF = (F x ,F y ,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 nondiagonal blocks as (F, K). For the system consisting of N nuclear spins, there are in total h = C N/2 N (binomial coefficient) diagonal blocks for the even N and h = C (N−1)/2 N 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 2 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ẑ (termed hereafter zero-quantum operators and denoted ZQ (F,K) ϕ,m ), a generalized measurement operator in a laboratory frame is written aŝ where the operators are defined up to a phase ϕ and a projection quantum number m such that |m| ≤ min(F, K).
When it comes to a diagonal block (F, F ), one can choose a measurement operator ZQ (F,F ) π/4,F , which is the same as the maximum-projection operator |F, F F, F |. The function of the surface representing the (F, F ) block can thus be expressed using Eq. (1): (F,F ) π/4,F rn = Tr ρRẑ →n |F, F F, F |R −1 z→n .
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); (F,F ) π/4,F rn 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 ZQ (F,K) ϕ,m such that ϕ = 0 and ϕ = π/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): with ϕ being equal to 0 or π/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.
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-bystep 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 (|1, 1 1, 1|)n measurements fully represent the related density operators (trivially, the (0, 0) block is represented by a sphere with a radius 1/4 and therefore is not shown). In Fig. 2c-d, we give an example of a density operatorρ =(1/4)1+(1/4)(Î 1z −Î 2z ) that is often encountered in various experiments involving two- Figure 2. (a-c) Visualizations using the measurement observable (|1, 1 1, 1|)n for the following density operators: (a)ρ = (1/4)1+(1/4)(Î 1z +Î 2z ), representing a state orientated alongẑ; , 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 ZQ (1,0) 0,0 n . spin systems (see the discussion below). Figure 2c shows the visualized surface based on the measurement (|1, 1 1, 1|)n, which appears isotropic. However, the state is intrinsically anisotropic along theẑ axis, which can be displayed by using a measurement operator ZQ (1,0) 0,0 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.

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 rn corresponding to surfaces constructed via Eqs. (5) and (6) contain as much information as the density matrix.
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),ρ =(1/4)1+Î 1zÎ2z , is shown in Fig. 3a; the density operator is fully represented by the surface when plotted using the measurement observable (|1, 1 1, 1|)n. Figure 3b shows the same density rotated by −(π/4) alongx.
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 surfaceplotted (|1, 1 1, 1|)n representing the state is 2-fold symmetric alongẑ. The related density matrix (written in Zeeman basis) only has nonzero coherences with | m| = 0 and | m| = 2. Figure 4b shows an example where two surfaces are needed to completely represent the state. The one plotted with the measurement operator (|1, 1 1, 1|)n (left) is 2fold symmetric aroundẑ, while that plotted with ZQ (1,0) π 2 ,0 n (middle) is rotationally symmetric (q-fold symmetric for an arbitrary positive integer q) aroundẑ. As a result, both surfaces possess a 2-fold symmetry aroundẑ, and the related density matrix only has nonzero coherences with | m| = 0 and | m| = 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 13 C-labeled formic acid with 1 H and 13 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, ν = (ν1 H + ν13 C )/2 (ν1 H and ν13 C are the 1 H and 13 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.
The observable operator describing a ZULF measurement is proportional to γ13 CÎ 1z + γ1 HÎ 2z (the operatorsÎ 1 andÎ 2 denote the 13 C and 1 H spins, respectively, and their projections are accordinglyÎ 1z andÎ 2z ). Further, one can decompose the observable operator into a symmetric  part (γ13 C + γ1 H )/2 (Î 1z +Î 2z ) and an antisymmetric part (γ13 C − γ1 H )/2 (Î 1z −Î 2z ). One may notice that symmetric and antisymmetric parts are proportional to the rank-1 component of the measurement operator |1, 1 1, 1| and the measurement operator ZQ (1,0) 0,0 , respectively. As a result, the ZULF spectra could be understood by checking the evolution of the surfaces plotted with (|1, 1 1, 1|)n and ZQ (1,0) 0,0 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ρ 0 =(1/4)1+ (1/10)(Î 1z + 4Î 2z ) (we assume high polarizations and set them proportional to gyromagnetic ratios using γ1 H /γ13 C ≈ 4) by plotting the surface with the operator (|1, 1 1, 1|)n. The resulting oriented surface precesses about the x axis with a period τ 1 = 1/ν. 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 ν, 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).
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,ρ 0 =(3/20)(Î 2z −Î 1z ), and symmetric,ρ 0 = (1/4)1+(1/4)(Î 1z +Î 2z ), parts of the initial density operator, respectively. One can see that the evolution of the asymmetric part of the density operator predominantly contributes to Figure 5. Visualizations of spin dynamics in an AX system ( 1 H-13 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 iŝ ρ 0 =(1/4)1+(1/10)(Î 1z +4Î 2z ) (Î 1 andÎ 2 denote 13 C and 1 H nuclei, respectively). Panels (c)-(d) show surfaces representing spin evolution in the ZULF experiment plotted with the measurements (c) (|1, 1 1, 1|)n over a timescale τ 1 = 1/ν = 2/(ν1 H + ν13 C ). Panel (d) ZQ (1,0) 0,0 n over a timescale τ 1 (inset) shows the evolution of the surfaces plotted in panel (d) over a shorter timescale τ 2 = 1/J . 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 ZQ (1,0) 0,0 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π J t) and the angle between the p orbital and the z axis being equal to (π + 2π νt) -the intersection of the surface with the z axis is found to have the following time dependence, cos(2π J t) cos(π + 2πνt) = −(1/2) [cos (2π (J + ν)t) + cos (2π(J − ν)t)], 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 π/2 radian out of phase in the orientation and −π/2 radian out of phase in the size oscillation. This means that the intersection contributing to the measurement has a time dependence cos(−π/2 + 2π J t) cos(3π/2 + 2πνt) = −(1/2) [cos (2π (J + ν)t) − cos (2π (J − ν)t)]. 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 ZQ (1,0) π/2,0 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.

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, = (δ 1 − δ 2 )γ B 0 , gives rise to singlet-triplet mixing. Such mixing is selective, and only |S → |T 0 transitions occur, while the other transitions |S → |T ± 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). 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 2 , respectively, and the longitudinal relaxation time of the two protons is given by T f 1 , when they are free in solution, and T b 1 , 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 → |T 0 transition and reduces the population difference between the two states. As the |T 0 state is populated, relaxation between triplet sublevels becomes more noticeable, and the states |T ± 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 0 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 ZQ (1,0) π 2 ,0 n , reveals extra information not covered by Fig. 7b; i.e., the singlet-triplet coherence,Î 1yÎ2x −Î 1xÎ2y , 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. 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 = √ 2/ . 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.

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 where ρ (F,K) λµ is the polarization moment and Y λµ (θ, φ) is the spherical harmonic. The defined radius can be complex. The amplitudes of (F,K) r(θ, φ) 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.
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 2 blocks: h = C N/2 N (binomial coefficient) when N is even and h = C (N −1)/2 N 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 + 1)/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 ν. 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 ν. While the surface also captures the full dynamics and directly reflects the two characteristic frequencies, J and ν, of the system, it is less obvious how to quantify the measured ZULF J spectra from it.

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 1 H-13 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 ultralowfield conditions.

Appendix A: Proof of rotational invariance ofÔn aroundn
Let us assume that a rotation about theẑ axis is applied to the measurement observableÔẑ. Since the rotation keeps the directionẑ unchanged, according to Eq. (2) the rotated observable can also be used for calculating rẑ. As a result, the two operators for measuring rẑ must be the same; otherwise, there would be ambiguity for the definition of rẑ.
Similarly, the observableÔn must be invariant under rotations aboutn (denotedRn) to avoid the ambiguity of rn. To see this, one can check the following equality: where we have applied Eq.
(2) and inserted the operator Rẑ →nR −1 z→n , which equals the identity on both sides. GivenR −1 z→nRnRẑ →n =Rẑ, one can show that which proves the conclusion.

B1 Definition of the spherical tensor basis operators
For each block (F, K), a basis of tensor operators is defined bŷ with |F − K| ≤ λ ≤ F + K and −λ ≤ µ ≤ λ and m K , m F , and µ being the projection quantum numbers of the Clebsch-Gordan coefficient Km K λµ|F m 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.

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

Hermitian adjoint
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. giving The sum can be evaluated as

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 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): [F z ,T λµ ] = µT λµ , Proof of the first equation is obvious. For the second part, we have We can apply Eq. (B4) and the recursion relationships between Clebsch-Gordan coefficients to find that 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: Given Eq. (B4) and the symmetry of the Clebsch-Gordan coefficients we have

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: Applying the orthogonality condition of the basis operators (Appendix B2), we have the so-called polarization moments (Auzinsh et al., 2010):

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 operatorRẑ →n asR(φ, θ, 0) (see Appendix B2 for a definition and assumen is pointing at (θ, φ)), the surfaces can be expressed using the function f There are two useful proprieties of the defined function f (F,K) m .

Complex conjugation
The first propriety (see the proof below) allows us to extract the polarization moments of the density operator from the function f Also, we note here that the second propriety guarantees that the surfaces plotted will always be real.
where we have employed the sphericalness of the basis operators (Appendix B2). The term F, m|T (F,K) λµ |K, m 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 Starting from Eq. (D2), by exchanging K and F , one can find that To find out the relationship between ρ By inserting and 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ẑ, and then all the elements ofρ (F,K) m F ,m 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 aŝ Following from the sphericalness (Appendix B2) of the basis operators, the polarization moments ofρ are Inserting ρ (F,K) λµ into Eq. (D2), one can derive the function f (F,K) m defined with respect to the rotated density operator ρ : On the other hand, one can apply the same rotation to surface Following from the rotational symmetry of spherical harmonics, one can realize f , 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., This requires all the elements ofρ (F,K) m F ,m K to be zero except those for which m F − m K = qN , where N is an integer.

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Ôẑ for AMPS is |F, F F, F |, which can be expanded over the spherical tensor operator basis (see Appendix B1) as |F, F F, F | = λ,µ c λ,µT where the coefficient c λ,µ can be calculated using Eq. (C2): c λ,µ = Tr |F, F F, F | T (F,K) The Clebsch-Gordan coefficient F F λµ|F F would only be nonzero when 0 ≤ λ ≤ 2F . When λ is within this range, one may find Inserting it into Eq. (G2), which is the same as the operator used for calculating the Husimi Q function (Koczor et al., 2020, Eq. 5).