# What Are the Relative Intensities of the Components of NMR Spectral Multiplets from Quadrupolar Nuclei in Uniformly Anisotropic Media?

Authors: Stuart J. Elliott; Philip W. Kuchel
Journal: Concepts in Magnetic Resonance Part A, Bridging Education and Research
 (2021)

---

## Abstract

Stretched hydrogels make uniformly anisotropic environments for quadrupolar nuclei such as2H, 23Na, and 133Cs. Such surroundings cause the partial alignment of nuclear spin bearing ions and molecules that is sufficiently pronounced to alter the nuclear magnetic resonance spectra of the guest species. In most cases, resonance splittings are directly related to the spin quantum number I. The relative intensities of the components of the resonance multiplets can be inferred from basic quantum mechanics.


---

## Body

## 1. Aim

The aim of this article is to undertake a step-by-step evolution of the formal mathematics that define the quantum mechanical description of the nuclear magnetic resonance (NMR) spectra of nuclei that have a quadrupole moment, when they are embedded within a molecularly anisotropic environment.


## 2. Motivation

This comes from viewing the well-defined septet obtained from133Cs+, with 133Cs NMR spectroscopy, when the cation is present in gelatin (and some other types of) gel, and the whole sample is held in a stretched state. Figure 1 shows such typical spectra.Figure 1
133Cs NMR (52.8 MHz) spectra (solid black lines) of 800 mM CsCl in 60% w/w gelatin gel, at 15°C, that was stretched by a factor of 2.0 (processed from 4 transients with 1 Hz spectral line broadening [1]). The splitting between adjacent peaks was 43.8 ± 0.1 Hz. (a) Spectral simulation (dashed magenta line) using a Markov chain Monte Carlo analysis and (b) the same experimental data as in (a), but the simulation (dashed red line) was based on an ensemble average of Gaussian distributed values of the peak splitting (reproduced with permission from the authors and publisher [1]).
(a)(b)A typical sample is made as follows: liquid gelatin (agarose or more exotic polymers) at 37–80°C containing the guest molecule or inorganic salt of interest (or even cells but not above 40°C) is drawn via an attached syringe into a silicone rubber tube. A Delrin (plastic) plug is inserted in the bottom end of the silicone tube which is then inserted into a bottomless thick-wall glass NMR tube. At the top end of the silicone tube, which projects beyond the glass tube, a thumb screw is attached. If the sample is to be compressed, then the silicone tube is stretched to whatever extent is prescribed, typically by a factor of 1.7 of the original length, and the thumb screw is tightened to hold this extent of stretching. The tube assembly is then placed in a refrigerator at 4°C for the gel to set, which it does once the temperature drops below ∼25°C. The sample is then compressed by releasing the thumb screw and allowing the silicone tube to relax fully or to be held at some intermediate extent of relaxation. For samples that are to be stretched, the tube assembly is placed in the refrigerator at 4°C without prior stretching, and stretching is carried out after the gel has set. The NMR spectra are recorded with the samples in the gel state, and this is typically at 20°C for gelatin gels [2–4].Why guest ions and molecules might be placed within gels and studied by NMR spectroscopy in the first place is a pertinent question. The motivation is to develop the growing use of reproducibly ordered media in analytical NMR spectroscopy [2, 5, 6]. The stretched/compressed gel approach is used in discriminating one enantiomer from another in chiral molecules [3, 7–10], identification and spectral separation of 13C-labelled isotopomers/isotopologues for the quantification of the relative proportions in a mixture [11], distinguishing the NMR signals from inside cells from those outside, comparing the ordering of the environments [1, 12], detecting heterogeneity by observing deviations from ideal spectral outcomes in gel-matrix systems [13], studying moiety mobility in small molecules [14], molecular structure refinement [15, 16], and most recently in exploring the effects on metabolic rate and membrane transport that are brought about by distorting cells, a phenomenon called molecular mechanosensation [4].It was during developing approaches for the latter experiments that the reproducible quadrupolar splittings in2H, 23Na, and 133Cs nuclei were discovered for stretched/compressed gelatin gels. The phenomenon became understood [17] and began to be exploited in studies of cells [1, 4].


## 3. Abstract Problem

At the abstract level, the problem is to compute the frequency distribution of the transitions between energy levels in a multilevel quantum system, when it relaxes from an excited state.


## 4. Communication

What constitutes a satisfactory explanation of a physical phenomenon depends on the technical and analytical background of the explainer and the recipient and on mutually understood language(s) and agreed symbolism. How the explanation links logically to the explanation of a related phenomenon often determines whether the explanation is considered to be good in the opinion of the recipient.


## 5. Basic Tenets

Here we define the various concepts in physics and then the elements of the mathematical machinery that are required to model these concepts, which in turn are required to solve the problem that is posed in the title of this article.


### 5.1. Quantization of Energy

This idea arose from the concept of quantization of light that emerged from Max Planck’s analysis of black body radiation and then Albert Einstein’s embracing of this idea as reality despite Planck’s initial protest that this was simply a conceptual aid [18].


### 5.2. Finite Number of Energy Levels

Another key idea is that for any atomic-scale system, its energy exists in discrete levels or states and that there are a finite number of these states. This deduction was based on optical spectroscopy of sun light in the early 1800s and formalized in 1913 into discrete energy levels in atoms by Niels Bohr [19]. It is not necessarily true that the number of possible quantum states, and hence number of transitions, is finite in a quantized atomic-scale system. This is exemplified by the hydrogen molecule in which the energy levels of the electrons are usually detected spectroscopically as finite in number, but in fact an electron can be extracted by various experimental means to an effectively infinite distance from the nucleus. Thus, it passes through an infinite number of energy levels. We will not develop this idea any further for nuclear spins except to say that the declaration of a finite number of energy levels fits with the sort of NMR spectroscopy we use in a terrestrial laboratory.


### 5.3. Processes of Excitation and Relaxation

The addition of energy to the physical system (excitation), or loss from it (relaxation), involves a redistribution of energy between the various states or a release of energy towards the lattice.


### 5.4. Rotation of a Vector by Rotation Matrices

This section is the run up to the so-called commutation relationships that exist between operators used in quantum mechanics, and as we will see, these apply to the matrix representations of operators as well. In their section on quantum mechanics, Hore et al. [20] state “Justifying (the commutation relationships) is not a simple business. One approach is to note that these forms (given in Part B of their book) give the correct commutation relationships.” So, it is this “not a simple business” that we will address here and develop in full.We start with classical vector analysis and consider how a matrix can be applied to a position vector (with its tail at the origin and tip at the position specified by the coordinate) to rotate the vector about the origin and hence shift its tip to another position. The three matrices that invoke rotations about thex-, y-, and z-axes of a Cartesian coordinate system are called the Euler rotation matrices. They have the respective forms denoted by [21](1)Rxθ=1000cosθ−sinθ0sinθcosθ,Ryθ=cosθ0sinθ010−sinθ0cosθ,Rzθ=cosθ−sinθ0sinθcosθ0001,where θ is the angle of rotation about the axis specified by the subscript. An illustration using Mathematica demonstrating the effect of these rotation matrices on a given vector is shown in Notebook 1.


### 5.5. Infinitesimal Rotations

Next consider the Euler matrices when there is an infinitesimal rotation, and for this step, we use the series expansions of the two trigonometric functions that appear above [21]:(2)sinε=ε−ε33!+ε55!−ε77!+⋯,cosε=1−ε22!+ε44!−ε66!+⋯.Hence, for small values ofε≪1, truncating the series after the cubic terms, which are even tinier, we obtain the following approximate forms of the Euler rotation matrices:(3)Rxε≈10001−ε22−ε0ε1−ε22,(4)Ryε≈1−ε220ε010−ε01−ε22,(5)Rzε≈1−ε22−ε0ε1−ε220001.Consider the single matrix that invokes an infinitesimal rotation about they-axis, followed by an infinitesimal rotation about the x-axis, noting that the operations are applied from the right; hence, we multiply the matrices as follows (where the dot denotes matrix multiplication):(6)Rxε⋅Ryε≈10001−ε22−ε0ε1−ε22.1−ε220ε010−ε01−ε22,=1−ε220εε21−ε22−ε−εε1−ε2.Now perform the rotations in reverse order:(7)Ryε⋅Rxε≈1−ε220ε010−ε01−ε22⋅10001−ε22−ε0ε1−ε22,=1−ε22ε2ε01−ε22−ε−εε1−ε2.Next, subtract the second product from the first to define a “quantification of the extent to which the order by which the operations are made actually matters mathematically.” Subtracting equation (7) from equation (6) gives(8)Rxε⋅Ryε−Ryε⋅Rxε≈0−ε20ε200000.The expression that arises from equation (8) is called the commutation relationship and will be used as a key concept throughout the remainder of the text (see later). Is equation (8) related to Rzε? Yes, except the latter has terms in ε, so let us replace these with their squares, which is valid since the square of an infinitesimal quantity remains infinitesimal. From equation (5), the expression for Rzθ with θ=ε2 gives(9)Rzε2≈1−ε42−ε20ε21−ε420001.And, with terms inε with powers greater than 2 sent to zero, we obtain(10)Rzε2≈1−ε20ε210001.This is decomposed into an interesting matrix plus the identity matrix:(11)Rzε2≈0−ε20ε200000+100010001,=0−ε20ε200000+1,where 1 denotes the identity matrix. Thus, the relationship in equation (8), which is called the commutation relationship, is given by(12)Rxε⋅Ryε−Ryε⋅Rxε=Rzε2−1.Also, note that the identity matrix can be viewed as a rotation of 0 radians about any axis, so we can write equation (12) as(13)Rxε⋅Ryε−Ryε⋅Rxε=Rzε2−Rany0.As can be seen in equation (13), the multiplication of the operators Rxε and Ryε, and their matrix representations, does not commute. We will return to this central expression later. In the meantime, let us develop some ideas about dynamical systems.


### 5.6. Infinitesimal Rotations of Angular Momentum

We will consider angular momentum first, since gaining ideas about it will allow us to transfer these to nuclear spin. The latter was a new concept back in the late 1920s; it is one that emerged from the classical idea of angular momentum, as we will see below.In general, we can represent a new operator,Uε (that in the present discussion can be thought of as a matrix, but this theory is quite general to operators), that arises from a tiny perturbation in its dependent variable, applied to the previous operator G:(14)Uε=I−iGε,where I denotes the identity operator (I is written in this skeletal form rather than a bold 1 to emphasize its operator nature rather than simply being a matrix as, say, in equation (12)), i is the imaginary unit (i = −1), and G is a general operator (considered to be Hermitian but what this is defined as and why it is important is considered below). In the present discussion, Uε is simply a perturbation of G, with the added insight that it is a complex function or operator, and hence generality is gained by multiplying it by i. We are thinking ahead here: we want our operator to be complex so that when we take an exponential of it, similar to eiθ and from Euler’s formula, we know this to be the sum of trigonometric functions, and hence it is periodic [21]: eiθ=cosθ+isinθ.In the case of angular momentum, we define the corresponding operatorJk (where k is the axis around which the rotation takes place; k = x, y or z) for an infinitesimal rotation of angle dϕ about the kth axis. And we replace Uε by a more evocative symbol D that means rotation (Drehung) in German. Thus, we can specify the z-axis and an infinitesimal angle of rotation by(15)Dzdϕ=I−iJzℏdϕ,where ℏ is simply the unit of energy that separates the finite number of different energy levels in our system. Note that we divide by ℏ (units J s, as expected for angular momentum) to normalize the term and to express the incremented angle in terms of the quantized energy.


### 5.7. Finite Rotations from Lots of Infinitesimal Ones

Now we apply a momentous mathematical insight. We use equation (15) to generalize to a finite angle of rotation, and for this, we use one of the standard definitions of “e” or “exp,” namely,(16)e=1+1nlimitn⟶∞n,or more generally:(17)e−x=1−xnlimitn⟶∞n.Thus, for a finite angle of rotationϕ, for which dϕ=ϕ/n, from equation (15), we obtain(18)Dzϕ=I−iJzℏϕnlimitn⟶∞n.So, from equation (17), we obtain(19)Dzϕ=exp−iJzℏϕ.Furthermore, we can expand the exponential to obtain(20)Dzϕ=I−iJzℏϕ+Jz/ℏϕ22!−⋯.And, by analogy with equations (13) and (20) becomes, on dropping higher degree terms:(21)Dzε2=I−iJzℏε2.By using our previous notation for incremental rotations,(22)Rzε2=I−iJzℏε2.Analogously, we write the expressions forRxε and Ryε using equation (20) coupled with truncating the series at the ε2 terms. So their products, according to equation (13), become(23)I−iJxεℏ−12!Jxεℏ2I−iJyεℏ−12!Jyεℏ2−I−iJyεℏ−12!Jyεℏ2I−iJxεℏ−12!Jxεℏ2=I−iJzℏε2−I.This is best shaken out inMathematica, noting that the order of the multiplications of the various Js must be preserved; otherwise, Mathematica will simply apply commutative multiplication and the left hand side will be found to be equal to 0 (see Notebook 1 for an implementation of this analysis in Mathematica). The tedious analysis gives the commutation relationship (see equation (8)):(24)JxJy−JyJx=iℏJz,or(25)Jx,Jy=iℏJz.And generally,(26)Ji,Jj=iℏεijkJk,where εijk is the so-called Levi-Civita permutation symbol which in two dimensions is(27)εij=+1,ifi,j=1,2−1,ifi,j=2,10,ifi=j.And, in three dimensions (with numbers 1, 2, and 3 arranged at the vertices of an equilateral triangle with 1 at the top, 2 at lower left, and 3 at lower right), using the cyclic diagram (↓1 ⟶ 2 ⟶ 3↑ meaning anticlockwise through the arrows) for +1 and (↑1 ⟵ 2 ⟵ 3↓ meaning clockwise through the arrows) for −1, we obtain(28)ɛijk=+1,ifi,j,k=1,2,3,2,3,1or3,1,2−1,ifi,j,k=1,3,2,2,1,3or3,2,10,ifi=j..Let us take an additional minor detour and really make sure we understand the parity calculation.Example 1.
Parity calculation and confirmation of equation (28).
Here is a permutation example inS5, or more simply put, a collection of 5 objects, where the top row in equation (29) is the argument of σ and the bottom row is the value of the function σ. Suppose(29)σ=1234545231.
The sign ofσ is the parity of the number of inversions when applying σ to the top row, i.e., the number of pairs (i, j) for which i < j such that σ (i) > σ (j). First, we have σ (1) = 4, σ (2) = 5 etc., so(30)σ1>σ3σ4σ5,σ2>σ3σ4σ5,σ3>σ5,σ4>σ5.
There are 8 inversions, and hencesgn (σ) = (−1)8 = 1. Thus, all the symbols used in equation (26), including the permutation symbol, are defined and ready for use in the theory we develop below to describe the spectrum (and those like it) in Figure 1.


### 5.8. Eigenvalues and Eigenstates of Angular Momentum

A general vector like angular momentum has an overall magnitude that is invariant under rotations, as discussed under equation (14). If it has components Jx, Jy, and Jz, then its magnitude is(31)J=JxJx+JyJy+JzJz,or, for ease of use in what follows, we can avoid the square root and instead define the operator as(32)J2=JxJx+JyJy+JzJz.It makes physical and mathematical sense that determining the magnitude of a vector (or operator) would not change if a rotation was to succeed or precede it, so it is no surprise thatJ2 commutes with Jk. Therefore,(33)J2,Jk=0,k=x,y,z.The consistency of this with the commutator expressions is seen by expanding equation (33) and using the commutation relationships [22] (which are relatively easy to prove) that are the same as the classical Poisson-bracket formulae of theoretical mechanics [23]. Equation (39) in particular will become very useful in what follows these:(34)A,A=0,(35)A,B=−B,A,(36)A,c=0,cis a scalar,(37)A+B,C=A,C+B,C,(38)A,BC=A,BC+BA,C,(39)AB,C=A,CB+AB,C,(40)A,B,C+B,C,A+C,A,B=0,Jacobi identity.Next, equation (33) can be written using equation (32) and the commutation relationships to confirm the physical reasoning given above:(41)JxJx+JyJy+JzJz,Jz=JxJx,Jz+JyJy,Jz+JzJz,Jz,=Jx,JzJx+JxJx,Jz+Jy,JzJy+JyJy,Jz+0,Jz,Jz=0,=−iℏJyJx−iℏJxJy+iℏJxJy+iℏJyJx,=0,where J2 serves as a conservation condition, and when it comes to the probability interpretation, i.e., expectation values of operators, we will anticipate that this will be normalized to 1 (but more on this later).


### 5.9. Simultaneous Eigenkets and the Raising and Lowering Operators

Eigenkets were first introduced by P. A. M. Dirac in order to put the whole of quantum mechanics on a more general footing. His novel notation to express typical quantum mechanical formulae uses bras and kets as a coordinate system for complex vector spaces [24]. The eigenkets contain all the information regarding the state of a quantum system and are represented by a column vector. This meaningful choice allows for an easy translation into whatever matrix representation is most useful for the problem at hand.It is reasonable to suppose that since the operatorsJx, Jy, and Jz do not commute with each other (equation (26)), we would expect that only one of them might share simultaneous eigenkets, i.e., be a commuting operator, with J2 (if at all; but suppose it does) and by convention this is chosen to be Jz. There should exist eigenkets and corresponding eigenvalues, and this is where we transition into quantum mechanics parlance and thinking, where we must postulate that the state of a system can be represented by state vectors and that operators operate on these to transform them. Also, the expectation of an operator, like angular momentum, relies on calculations with these eigenstates and eigenvectors [20, 22]. For the moment, consider(42)J2a,b=aa,b,Jza,b=ba,b,where the eigenket is labelled (given a name) based on whichever eigenvalues it yields when operated upon by J2 (namely, a) or Jz (namely, b). This convention is violated by many authors and we join them later (see below), but it is useful here to be clear about what a and b mean inside the bra-ket; they are simply labels.Thinking ahead, we want to develop the mathematics so that it models the physical fact that there are a finite number of energy levels and hence eigenvectors and that the magnitude (squared also) of the vectors and corresponding operators is fixed.We introduce a composite operator and explore its properties:(43)J±=Jx±iJy.There is a surprising and most useful property, from the point of view of the subsequent theory, that follows some algebraic manipulations. Thus, the following commutation relationships can be shown to hold:(44)J+,J−=2ℏJz,(45)Jz,J±=±ℏJ±.And from equation (33),(46)J2,J±=0.But things get interesting from the point of view of the physical meaning of theJ± operator. To expose this trickiness, consider the following equation that describes the sequential operation of Jz on J±, which is operating directly on the eigenket a,b. Thus, expanding equation (45) gives(47)JzJ±−J±Jza,b=±ℏJ±a,b,(48)JzJ±a,b=J±Jza,b±ℏJ±a,b,(49)JzJ±a,b=J±ba,b±ℏJ±a,b,(50)JzJ±a,b=b±hJ±a,b.Hence, ifJz acts on the eigenket that has been transformed by J±, then its eigenvalue is its original eigenvalue plus ℏ. Because of this raising +ℏ or lowering −ℏ of the energy eigenvalue, the operators J+ and J− are called the “raising” and “lowering” operators.SinceJ± and J2 commute, equation (46) specifies the following:(51)J2J±a,b=J±J2a,b,=J±aa,b.So,J±a,b is an eigenket of J2 and Jz. Therefore, using the aforementioned labelling convention for eigenvectors, we can write J±a,b=ca,b±ℏ, and we leave it to later to sort out what form c takes.


### 5.10. The Range of Eigenvalues ofJz

We seek a general expression for the eigenvalues of the eigenkets of the angular momentum operator. We can applyJ+ to the eigenkets but eventually a maximum eigenvalue is reached; call this maximum value nℏ. Of course, applying J2 to the eigenkets leaves the eigenvalue unchanged, as in equation (51). In other words, there is an upper limit to b for a given value of a in a,b.In fact, we shall prove thata≥b2 (meaning that the magnitude of the z-component of a vector cannot exceed that of the vector as a whole), so to begin the analysis, note that [22](52)J2−Jz2=Jx2+Jy2.But,(53)Jx=12J++J−,Jy=−i2J+−J−.Therefore,(54)J2−Jz2=12J+2+J−2.Physically, there is a maximum eigenvalue beyond which the raising operator cannot raise the eigenket, so let us call this eigenvaluebmax; in other words, J+a,bmax=0.Also, it makes physical sense that(55)J−J+a,bmax=0.But,(56)J−J+=Jx2+Jy2−iJyJx−JyJx,=Jx2+Jy2−iJy,Jx.From equation (52),(57)=J2−Jz2−iJy,Jx.And now we use the much vaunted commutator that we took so long to derive up to equation (25):(58)=J2−Jz2−ℏJz.When applied to the maximum eigenvalue eigenket, we obtain(59)J2−Jz2−ℏJza,bmax=a2−bmax2−ℏbmaxa,bmax=0.Sincea,bmax is not a null eigenket, i.e., it exists in general, for the zero to hold in general, then(60)a=bmaxbmax+ℏ.Sincebmax is the largest eigenvalue and given that ℏ is positive, then in general, a≥b2.Also, physically, there is a lower value of the energy level below which the operators cannot move the system. This needs to be built into the model of the system. We call this minimum eigenvaluebmin such that J−a,bmin=0. By analogy with equation (56),(61)J+J−=Jx2+Jy2+iJyJx−JyJx,=J2−Jz2+ℏJz,from which, using the counterpart of equation (59), we obtain(62)a=bminbmin−ℏ.Comparing equations (60) and (62), we see that(63)bmin=−bmax.In other words,(64)−bmax<b<bmax.By applyingJ+ successively to a,bmin a finite number of times, say n times, we have(65)bmax=bmin+nℏ.Thus,(66)bmax=nℏ2.It is conventional to work in units ofbmax/ℏ=jso thatj=n/2. Thus, the maximum value of the Jz eigenvalue is jℏ, where, as we will see later, j can be either an integer or half-integer. Furthermore, since bmax=jℏ, equation (60) can be written as(67)a=jℏjℏ+ℏ=ℏ2jj+1.As for the eigenvalue ofJ2, we can define b=mℏ, and by sorting through the algebra, it is readily shown that the allowed values of m are(68)m=−j,−j−1,−j−2,…,j−2,j−1,j,giving a total of2j+1states.Overall, instead of using the eigenvector designationa,b, we write the new form of the eigenvalues inside the ket not as ℏ2jj+1,mℏ but simply as j,m. Nevertheless, we must remember that(69)J2j,m=ℏ2jj+1j,m,Jzj,m=mℏj,m.However, we still have to find a value forc in the equation J±a,b=ca,b±ℏ.


## 5.1. Quantization of Energy

This idea arose from the concept of quantization of light that emerged from Max Planck’s analysis of black body radiation and then Albert Einstein’s embracing of this idea as reality despite Planck’s initial protest that this was simply a conceptual aid [18].


## 5.2. Finite Number of Energy Levels

Another key idea is that for any atomic-scale system, its energy exists in discrete levels or states and that there are a finite number of these states. This deduction was based on optical spectroscopy of sun light in the early 1800s and formalized in 1913 into discrete energy levels in atoms by Niels Bohr [19]. It is not necessarily true that the number of possible quantum states, and hence number of transitions, is finite in a quantized atomic-scale system. This is exemplified by the hydrogen molecule in which the energy levels of the electrons are usually detected spectroscopically as finite in number, but in fact an electron can be extracted by various experimental means to an effectively infinite distance from the nucleus. Thus, it passes through an infinite number of energy levels. We will not develop this idea any further for nuclear spins except to say that the declaration of a finite number of energy levels fits with the sort of NMR spectroscopy we use in a terrestrial laboratory.


## 5.3. Processes of Excitation and Relaxation

The addition of energy to the physical system (excitation), or loss from it (relaxation), involves a redistribution of energy between the various states or a release of energy towards the lattice.


## 5.4. Rotation of a Vector by Rotation Matrices

This section is the run up to the so-called commutation relationships that exist between operators used in quantum mechanics, and as we will see, these apply to the matrix representations of operators as well. In their section on quantum mechanics, Hore et al. [20] state “Justifying (the commutation relationships) is not a simple business. One approach is to note that these forms (given in Part B of their book) give the correct commutation relationships.” So, it is this “not a simple business” that we will address here and develop in full.We start with classical vector analysis and consider how a matrix can be applied to a position vector (with its tail at the origin and tip at the position specified by the coordinate) to rotate the vector about the origin and hence shift its tip to another position. The three matrices that invoke rotations about thex-, y-, and z-axes of a Cartesian coordinate system are called the Euler rotation matrices. They have the respective forms denoted by [21](1)Rxθ=1000cosθ−sinθ0sinθcosθ,Ryθ=cosθ0sinθ010−sinθ0cosθ,Rzθ=cosθ−sinθ0sinθcosθ0001,where θ is the angle of rotation about the axis specified by the subscript. An illustration using Mathematica demonstrating the effect of these rotation matrices on a given vector is shown in Notebook 1.


## 5.5. Infinitesimal Rotations

Next consider the Euler matrices when there is an infinitesimal rotation, and for this step, we use the series expansions of the two trigonometric functions that appear above [21]:(2)sinε=ε−ε33!+ε55!−ε77!+⋯,cosε=1−ε22!+ε44!−ε66!+⋯.Hence, for small values ofε≪1, truncating the series after the cubic terms, which are even tinier, we obtain the following approximate forms of the Euler rotation matrices:(3)Rxε≈10001−ε22−ε0ε1−ε22,(4)Ryε≈1−ε220ε010−ε01−ε22,(5)Rzε≈1−ε22−ε0ε1−ε220001.Consider the single matrix that invokes an infinitesimal rotation about they-axis, followed by an infinitesimal rotation about the x-axis, noting that the operations are applied from the right; hence, we multiply the matrices as follows (where the dot denotes matrix multiplication):(6)Rxε⋅Ryε≈10001−ε22−ε0ε1−ε22.1−ε220ε010−ε01−ε22,=1−ε220εε21−ε22−ε−εε1−ε2.Now perform the rotations in reverse order:(7)Ryε⋅Rxε≈1−ε220ε010−ε01−ε22⋅10001−ε22−ε0ε1−ε22,=1−ε22ε2ε01−ε22−ε−εε1−ε2.Next, subtract the second product from the first to define a “quantification of the extent to which the order by which the operations are made actually matters mathematically.” Subtracting equation (7) from equation (6) gives(8)Rxε⋅Ryε−Ryε⋅Rxε≈0−ε20ε200000.The expression that arises from equation (8) is called the commutation relationship and will be used as a key concept throughout the remainder of the text (see later). Is equation (8) related to Rzε? Yes, except the latter has terms in ε, so let us replace these with their squares, which is valid since the square of an infinitesimal quantity remains infinitesimal. From equation (5), the expression for Rzθ with θ=ε2 gives(9)Rzε2≈1−ε42−ε20ε21−ε420001.And, with terms inε with powers greater than 2 sent to zero, we obtain(10)Rzε2≈1−ε20ε210001.This is decomposed into an interesting matrix plus the identity matrix:(11)Rzε2≈0−ε20ε200000+100010001,=0−ε20ε200000+1,where 1 denotes the identity matrix. Thus, the relationship in equation (8), which is called the commutation relationship, is given by(12)Rxε⋅Ryε−Ryε⋅Rxε=Rzε2−1.Also, note that the identity matrix can be viewed as a rotation of 0 radians about any axis, so we can write equation (12) as(13)Rxε⋅Ryε−Ryε⋅Rxε=Rzε2−Rany0.As can be seen in equation (13), the multiplication of the operators Rxε and Ryε, and their matrix representations, does not commute. We will return to this central expression later. In the meantime, let us develop some ideas about dynamical systems.


## 5.6. Infinitesimal Rotations of Angular Momentum

We will consider angular momentum first, since gaining ideas about it will allow us to transfer these to nuclear spin. The latter was a new concept back in the late 1920s; it is one that emerged from the classical idea of angular momentum, as we will see below.In general, we can represent a new operator,Uε (that in the present discussion can be thought of as a matrix, but this theory is quite general to operators), that arises from a tiny perturbation in its dependent variable, applied to the previous operator G:(14)Uε=I−iGε,where I denotes the identity operator (I is written in this skeletal form rather than a bold 1 to emphasize its operator nature rather than simply being a matrix as, say, in equation (12)), i is the imaginary unit (i = −1), and G is a general operator (considered to be Hermitian but what this is defined as and why it is important is considered below). In the present discussion, Uε is simply a perturbation of G, with the added insight that it is a complex function or operator, and hence generality is gained by multiplying it by i. We are thinking ahead here: we want our operator to be complex so that when we take an exponential of it, similar to eiθ and from Euler’s formula, we know this to be the sum of trigonometric functions, and hence it is periodic [21]: eiθ=cosθ+isinθ.In the case of angular momentum, we define the corresponding operatorJk (where k is the axis around which the rotation takes place; k = x, y or z) for an infinitesimal rotation of angle dϕ about the kth axis. And we replace Uε by a more evocative symbol D that means rotation (Drehung) in German. Thus, we can specify the z-axis and an infinitesimal angle of rotation by(15)Dzdϕ=I−iJzℏdϕ,where ℏ is simply the unit of energy that separates the finite number of different energy levels in our system. Note that we divide by ℏ (units J s, as expected for angular momentum) to normalize the term and to express the incremented angle in terms of the quantized energy.


## 5.7. Finite Rotations from Lots of Infinitesimal Ones

Now we apply a momentous mathematical insight. We use equation (15) to generalize to a finite angle of rotation, and for this, we use one of the standard definitions of “e” or “exp,” namely,(16)e=1+1nlimitn⟶∞n,or more generally:(17)e−x=1−xnlimitn⟶∞n.Thus, for a finite angle of rotationϕ, for which dϕ=ϕ/n, from equation (15), we obtain(18)Dzϕ=I−iJzℏϕnlimitn⟶∞n.So, from equation (17), we obtain(19)Dzϕ=exp−iJzℏϕ.Furthermore, we can expand the exponential to obtain(20)Dzϕ=I−iJzℏϕ+Jz/ℏϕ22!−⋯.And, by analogy with equations (13) and (20) becomes, on dropping higher degree terms:(21)Dzε2=I−iJzℏε2.By using our previous notation for incremental rotations,(22)Rzε2=I−iJzℏε2.Analogously, we write the expressions forRxε and Ryε using equation (20) coupled with truncating the series at the ε2 terms. So their products, according to equation (13), become(23)I−iJxεℏ−12!Jxεℏ2I−iJyεℏ−12!Jyεℏ2−I−iJyεℏ−12!Jyεℏ2I−iJxεℏ−12!Jxεℏ2=I−iJzℏε2−I.This is best shaken out inMathematica, noting that the order of the multiplications of the various Js must be preserved; otherwise, Mathematica will simply apply commutative multiplication and the left hand side will be found to be equal to 0 (see Notebook 1 for an implementation of this analysis in Mathematica). The tedious analysis gives the commutation relationship (see equation (8)):(24)JxJy−JyJx=iℏJz,or(25)Jx,Jy=iℏJz.And generally,(26)Ji,Jj=iℏεijkJk,where εijk is the so-called Levi-Civita permutation symbol which in two dimensions is(27)εij=+1,ifi,j=1,2−1,ifi,j=2,10,ifi=j.And, in three dimensions (with numbers 1, 2, and 3 arranged at the vertices of an equilateral triangle with 1 at the top, 2 at lower left, and 3 at lower right), using the cyclic diagram (↓1 ⟶ 2 ⟶ 3↑ meaning anticlockwise through the arrows) for +1 and (↑1 ⟵ 2 ⟵ 3↓ meaning clockwise through the arrows) for −1, we obtain(28)ɛijk=+1,ifi,j,k=1,2,3,2,3,1or3,1,2−1,ifi,j,k=1,3,2,2,1,3or3,2,10,ifi=j..Let us take an additional minor detour and really make sure we understand the parity calculation.Example 1.
Parity calculation and confirmation of equation (28).
Here is a permutation example inS5, or more simply put, a collection of 5 objects, where the top row in equation (29) is the argument of σ and the bottom row is the value of the function σ. Suppose(29)σ=1234545231.
The sign ofσ is the parity of the number of inversions when applying σ to the top row, i.e., the number of pairs (i, j) for which i < j such that σ (i) > σ (j). First, we have σ (1) = 4, σ (2) = 5 etc., so(30)σ1>σ3σ4σ5,σ2>σ3σ4σ5,σ3>σ5,σ4>σ5.
There are 8 inversions, and hencesgn (σ) = (−1)8 = 1. Thus, all the symbols used in equation (26), including the permutation symbol, are defined and ready for use in the theory we develop below to describe the spectrum (and those like it) in Figure 1.


## 5.8. Eigenvalues and Eigenstates of Angular Momentum

A general vector like angular momentum has an overall magnitude that is invariant under rotations, as discussed under equation (14). If it has components Jx, Jy, and Jz, then its magnitude is(31)J=JxJx+JyJy+JzJz,or, for ease of use in what follows, we can avoid the square root and instead define the operator as(32)J2=JxJx+JyJy+JzJz.It makes physical and mathematical sense that determining the magnitude of a vector (or operator) would not change if a rotation was to succeed or precede it, so it is no surprise thatJ2 commutes with Jk. Therefore,(33)J2,Jk=0,k=x,y,z.The consistency of this with the commutator expressions is seen by expanding equation (33) and using the commutation relationships [22] (which are relatively easy to prove) that are the same as the classical Poisson-bracket formulae of theoretical mechanics [23]. Equation (39) in particular will become very useful in what follows these:(34)A,A=0,(35)A,B=−B,A,(36)A,c=0,cis a scalar,(37)A+B,C=A,C+B,C,(38)A,BC=A,BC+BA,C,(39)AB,C=A,CB+AB,C,(40)A,B,C+B,C,A+C,A,B=0,Jacobi identity.Next, equation (33) can be written using equation (32) and the commutation relationships to confirm the physical reasoning given above:(41)JxJx+JyJy+JzJz,Jz=JxJx,Jz+JyJy,Jz+JzJz,Jz,=Jx,JzJx+JxJx,Jz+Jy,JzJy+JyJy,Jz+0,Jz,Jz=0,=−iℏJyJx−iℏJxJy+iℏJxJy+iℏJyJx,=0,where J2 serves as a conservation condition, and when it comes to the probability interpretation, i.e., expectation values of operators, we will anticipate that this will be normalized to 1 (but more on this later).


## 5.9. Simultaneous Eigenkets and the Raising and Lowering Operators

Eigenkets were first introduced by P. A. M. Dirac in order to put the whole of quantum mechanics on a more general footing. His novel notation to express typical quantum mechanical formulae uses bras and kets as a coordinate system for complex vector spaces [24]. The eigenkets contain all the information regarding the state of a quantum system and are represented by a column vector. This meaningful choice allows for an easy translation into whatever matrix representation is most useful for the problem at hand.It is reasonable to suppose that since the operatorsJx, Jy, and Jz do not commute with each other (equation (26)), we would expect that only one of them might share simultaneous eigenkets, i.e., be a commuting operator, with J2 (if at all; but suppose it does) and by convention this is chosen to be Jz. There should exist eigenkets and corresponding eigenvalues, and this is where we transition into quantum mechanics parlance and thinking, where we must postulate that the state of a system can be represented by state vectors and that operators operate on these to transform them. Also, the expectation of an operator, like angular momentum, relies on calculations with these eigenstates and eigenvectors [20, 22]. For the moment, consider(42)J2a,b=aa,b,Jza,b=ba,b,where the eigenket is labelled (given a name) based on whichever eigenvalues it yields when operated upon by J2 (namely, a) or Jz (namely, b). This convention is violated by many authors and we join them later (see below), but it is useful here to be clear about what a and b mean inside the bra-ket; they are simply labels.Thinking ahead, we want to develop the mathematics so that it models the physical fact that there are a finite number of energy levels and hence eigenvectors and that the magnitude (squared also) of the vectors and corresponding operators is fixed.We introduce a composite operator and explore its properties:(43)J±=Jx±iJy.There is a surprising and most useful property, from the point of view of the subsequent theory, that follows some algebraic manipulations. Thus, the following commutation relationships can be shown to hold:(44)J+,J−=2ℏJz,(45)Jz,J±=±ℏJ±.And from equation (33),(46)J2,J±=0.But things get interesting from the point of view of the physical meaning of theJ± operator. To expose this trickiness, consider the following equation that describes the sequential operation of Jz on J±, which is operating directly on the eigenket a,b. Thus, expanding equation (45) gives(47)JzJ±−J±Jza,b=±ℏJ±a,b,(48)JzJ±a,b=J±Jza,b±ℏJ±a,b,(49)JzJ±a,b=J±ba,b±ℏJ±a,b,(50)JzJ±a,b=b±hJ±a,b.Hence, ifJz acts on the eigenket that has been transformed by J±, then its eigenvalue is its original eigenvalue plus ℏ. Because of this raising +ℏ or lowering −ℏ of the energy eigenvalue, the operators J+ and J− are called the “raising” and “lowering” operators.SinceJ± and J2 commute, equation (46) specifies the following:(51)J2J±a,b=J±J2a,b,=J±aa,b.So,J±a,b is an eigenket of J2 and Jz. Therefore, using the aforementioned labelling convention for eigenvectors, we can write J±a,b=ca,b±ℏ, and we leave it to later to sort out what form c takes.


## 5.10. The Range of Eigenvalues ofJz

We seek a general expression for the eigenvalues of the eigenkets of the angular momentum operator. We can applyJ+ to the eigenkets but eventually a maximum eigenvalue is reached; call this maximum value nℏ. Of course, applying J2 to the eigenkets leaves the eigenvalue unchanged, as in equation (51). In other words, there is an upper limit to b for a given value of a in a,b.In fact, we shall prove thata≥b2 (meaning that the magnitude of the z-component of a vector cannot exceed that of the vector as a whole), so to begin the analysis, note that [22](52)J2−Jz2=Jx2+Jy2.But,(53)Jx=12J++J−,Jy=−i2J+−J−.Therefore,(54)J2−Jz2=12J+2+J−2.Physically, there is a maximum eigenvalue beyond which the raising operator cannot raise the eigenket, so let us call this eigenvaluebmax; in other words, J+a,bmax=0.Also, it makes physical sense that(55)J−J+a,bmax=0.But,(56)J−J+=Jx2+Jy2−iJyJx−JyJx,=Jx2+Jy2−iJy,Jx.From equation (52),(57)=J2−Jz2−iJy,Jx.And now we use the much vaunted commutator that we took so long to derive up to equation (25):(58)=J2−Jz2−ℏJz.When applied to the maximum eigenvalue eigenket, we obtain(59)J2−Jz2−ℏJza,bmax=a2−bmax2−ℏbmaxa,bmax=0.Sincea,bmax is not a null eigenket, i.e., it exists in general, for the zero to hold in general, then(60)a=bmaxbmax+ℏ.Sincebmax is the largest eigenvalue and given that ℏ is positive, then in general, a≥b2.Also, physically, there is a lower value of the energy level below which the operators cannot move the system. This needs to be built into the model of the system. We call this minimum eigenvaluebmin such that J−a,bmin=0. By analogy with equation (56),(61)J+J−=Jx2+Jy2+iJyJx−JyJx,=J2−Jz2+ℏJz,from which, using the counterpart of equation (59), we obtain(62)a=bminbmin−ℏ.Comparing equations (60) and (62), we see that(63)bmin=−bmax.In other words,(64)−bmax<b<bmax.By applyingJ+ successively to a,bmin a finite number of times, say n times, we have(65)bmax=bmin+nℏ.Thus,(66)bmax=nℏ2.It is conventional to work in units ofbmax/ℏ=jso thatj=n/2. Thus, the maximum value of the Jz eigenvalue is jℏ, where, as we will see later, j can be either an integer or half-integer. Furthermore, since bmax=jℏ, equation (60) can be written as(67)a=jℏjℏ+ℏ=ℏ2jj+1.As for the eigenvalue ofJ2, we can define b=mℏ, and by sorting through the algebra, it is readily shown that the allowed values of m are(68)m=−j,−j−1,−j−2,…,j−2,j−1,j,giving a total of2j+1states.Overall, instead of using the eigenvector designationa,b, we write the new form of the eigenvalues inside the ket not as ℏ2jj+1,mℏ but simply as j,m. Nevertheless, we must remember that(69)J2j,m=ℏ2jj+1j,m,Jzj,m=mℏj,m.However, we still have to find a value forc in the equation J±a,b=ca,b±ℏ.


## 6. Elements of the Matrix Representation of the Angular Momentum Operator

Hore et al. [20] and Sakuri and Napolitano [22] give instructions on how to represent an operator as a matrix, composed using a set of eigenkets (and we expand on this in an Aside starting just before equation (80)). We take up the story from Sakuri and Napolitano [22] and suppose that the eigenkets of the basis set are normalized (sum their eigenvalues and divide each separate eigenvalue by this sum) in which case(70)j′,m′J2j,m=jj+1ℏ2δj′jδm′m,where δj′j and δm′m are Dirac delta functions, and(71)j′,m′Jzj,m=mℏ2δj′jδm′m.To obtain the matrix elements ofJ±, we use equation (58) in combination with equations (70) and (71):(72)j,mJ−J+j,m=j,mJ2−Jz2−ℏJzj,m,=ℏ2jj+1−m2−m.It follows from the comments under equation (51) that there exists a value of cjm+ which satisfies(73)J+j,m=cjm+j,m+1.So, it follows from equation (72) that(74)cjm+2=ℏ2jj+1−m2−m.And substituting equation (74) into equation (73), we obtain(75)J+j,m=ℏj−mj+m+1j,m+1.Similarly,(76)J−j,m=ℏj+mj−m+1j,m−1.Finally, the matrix elements of the matrix representation ofJ± are given by(77)j′,m′J±j,m=j∓mj±m+1ℏδj′jδm′m±1,i.e., the value of c which we sought above is given by j∓mj±m+1. Finally, this expression is the means whereby we obtain all the elements of the matrix representations of the nuclear spin operators. Therefore, it is important to consider what we mean by nuclear spin.


## 7. Spin

The theoretical developments for nuclear spin are exactly the same as for the more conventional angular momentum given above. In fact, we can simply substituteI2, Ix, Iy, and Iz wherever J2, Jx, Jy, and Jz appear in the equations.Specifically, what remains to be done in the case of the main objective of this article is to use the spin version of equation (77) to write the matrix representation of I±, each term of which is(78)I′,m′I±I,m=I∓mI±m+1ℏδI′Iδm′m±1.After making the substitutionm = −I, −I + 1, −I + 2, …, I − 2, I − 1, I (equation (68)), it becomes clear that the matrix of the lowering operator (equation (78)) has the following form:(79)I−=000⋯00002I0000000022I−1000000006I−100000⋮00⋱000⋮00006I−10000000022I−100000⋯002I0,see Notebook 2 for the sequence of events used to generate equation (79).Equation (79) highlights the complicated dependence of the matrix representation for the lowering operator I− on the value of the nuclear spin angular momentum quantum number I. The cases for I = 1, 3/2, and 7/2 are very much relevant to the problem of determining spectral multiplet patterns in the NMR spectroscopy of quadrupolar nuclei embedded within stretched hydrogel media. These examples are described in detail in the following sections. There we will discover how the matrix form of the lowering operator in equation (79) feeds into the multiplet structure of the acquired NMR spectrum for a nucleus with spin angular momentum quantum number I.The pertinent questions at this point in our exploration of NMR quantum mechanics are as follows: (i) What state do we desire the spin system to possess just before signal detection? (ii) And how is this target state generated in an NMR experiment? To address these important practical questions, we need another subtle mathematical object.


## 8. General Matrix Representation of an Operator

Aside: in all the sections below, we apply a caret to those symbols that denote operators. This was not necessary above because essentially all the dependent variables were operators and no confusion should arise since we know that operator algebra applies there. However, in the developments below, we will ultimately want to introduce the notion of superoperators, and for them, we need an identifier to distinguish them from simple operators; we use a double caret for this.The general matrix representation of a quantum mechanical operator is clearly described by Hore et al. [20]: if the operator A^ acts on the eigenket ψ with eigenvalue a, then(80)A^ψ=aψ.IfA^ corresponds to an experimentally observable property of the system, such as energy or angular momentum, then the mean value is defined by the expectation value, which is simply the scalar product of ψ with A^ψ that is denoted by(81)A^=ψA^ψ,where ψ=ψ† and † indicates the conjugate transpose. In general, the basis states of a system are not eigenstates of A^, e.g., for a two-state system with basis states 1 and 2, where ψ=c11+c22, we find that(82)A^=1c1∗+2c2∗A^c11+c22,=c1∗c11A^1+c1∗c21A^2+c2∗c12A^1+c2∗c22A^2.And this is the key insight: we can consider each termiA^j in equation (82) as an element Aij in the matrix representation of A^, which we denote in bold font style:(83)A=A11A12A21A22=1A^11A^22A^12A^2.Hence, equation (82) can be written as(84)A^=c1∗c1A11+c1∗c2A12+c2∗c1A21+c2∗c2A22.With this basic idea in mind, we can generalize the approach to higher dimensions and make compact matrix representations of multiple spin systems.


## 9. Spin Density Operator

In general, the state of a single spin can be described by a superposition state, i.e., the quantum state of a spinψ can be expanded according to a linear combination of suitable basis kets denoted by i:(85)ψ=∑icii,where ci is a complex coefficient (different from c discussed in equations (51) and (69)).As shown in equations (80)–(84), the expectation value of an operator depends on the products of coefficients. As a result, it becomes convenient to define a “spin density operator” ρ^ which can be described easily in terms of the elements of a matrix, e.g., between states i and j:(86)ρij=iρ^j=cicj∗.This facilitates the construction of the matrix form ofρ^, which is expressed in terms of complex coefficients:(87)ρ^=cicjcicj,=cici∗cicj∗cjci∗cjcj∗.As such, equation (84) can be reformulated in a more compact manner:(88)A^=ρ11A11+ρ12A12+ρ21A21+ρ22A22.Equation (87) shows how the spin density matrix can be constructed from underlying coefficients. Therefore, by using the analogies between basis vectors and eigenkets, ρ^ can be reformulated as(89)ρ^=ψψ,=∑i,jcicj∗ij.This relationship suggests that for a specific state of a spinψ, the quantity ψρ^ψ provides the probability of finding the spin in that state. Specifically, a foundational concept in quantum mechanics is that the expectation value of an operator A^ acting on ψ may be extracted as follows (in agreement with equation (81)):(90)A^=∑i,jcicj∗iA^j.By using the substitution in equation (86), and similar to equation (83):(91)Aij=iA^j.Equation (90) simplifies to(92)A^=∑ijρijAij,=Trρ^A^,where Tr is the trace operation (sum of the diagonal elements) which acts over a matrix [20, 25].It is worth noting that the total magnetization for aspin ensemble is the sum over all nuclear magnetic moments of the individual spins. This summation would be exceedingly tedious to do considering that an ensemble (in an NMR sample tube) typically has ∼1022 spins, say, in a solution of 0.1 mol L−1spin-bearing solute. However, given that the mathematics developed above can be used to encapsulate the dynamics of each member within the whole nuclear ensemble, then equation (89) can be ensemble averaged, and thus ρ^ reports on the overall spin state of the ensemble without explicitly referring to individual spins. How this idea is turned into mathematics is very subtle. It is superbly described in Hore et al. [20], Levitt [25], and Keeler [26].The physical interpretation of the spin density operatorρ^ is remarkably simple. It posits that the result of any repeated observations of identically prepared systems can be determined by using two properties which are related to quantum mechanical operators: (i) the state of the entire spin system (given by ρ^); and (ii) the property under measurement (given by an operator, e.g., I^−). The result is independent of the number of spins in the system, which provides a great simplification compared to specifying the states of ∼1022 spins. ρ^ is therefore related to a quantum probability density, but it is unfortunate that the name bears no physical relationship to the number density of nuclear spins within a sample.It is now reasonable to ask “How do we get away with making this bold generalization and get something useful out of the analysis?” One illustration is that equation (92) can be used to calculate the z-component of a sample’s magnetization M^z:(93)M^z=Trρ^M^z,=Nγℏ×TrI^zρ^,where N is the number of spins in the sample and γ is their magnetogyric ratio. For a noninteracting ensemble of nuclei at thermal equilibrium, the trace operation yields the expected result for the sample’s z-magnetization component:(94)M^z=Nγ2ℏ2B04kBT,where B0 is the static magnetic field strength and T is the temperature. The significance of the expectation value of the lowering operator in determining the spectral outcome of an NMR experiment is discussed in the next sections.


## 10. Liouville–von Neumann Equation

The time evolution of the state of a spin systemψt is represented by an equation of motion, known as the Schrödinger equation:(95)ddtψt=−iℏH^tψt,where H^t is the so-called Hamiltonian or total energy operator. The Hamiltonian specifies how the ensemble spin state ψt evolves in time.The Liouville–von Neumann equation is directly derived from the Schrödinger equation using the definition ofρ^t (see equation (89)) and differentiating by parts:(96)ddtρ^t=ddtψtψt,=ddtψtψt+ψtddtψt.A more intuitive derivation of the Liouville–von Neumann equation uses the outer product withψt in the Schrödinger equation (Equation (95)), thus identifying a commutator relationship:(97)ddtρ^t=−iℏH^tψtψt+iℏψtψtH^t,=−iℏH^t,ρ^t.Hence, the equation of motion of the density operatorρ^t is the Liouville–von Neumann equation. The solution of it is given by(98)ρ^t=U^t,0ρ^0U^t,0†,where † indicates the conjugate transpose and U^t,0 is a time-propagation operator:(99)U^t,0=e−iℋ^t/ℏ.The veracity of equation (98) as the solution of equation (97) is shown simply by differentiating it. Equation (98) neatly allows for the prediction of the state of the spin system at a later time given the value at t = 0 and thus describes the evolution of ρ^t under the given Hamiltonian operator H^.The commutation relationship at the heart of equation (97) can be used to infer well-known properties of spin systems in magnetic resonance experiments, e.g., if the state of the spin system commutes with the Hamiltonian, then the system does not evolve (in the absence of nuclear relaxation). Practical examples describing the functionality of the Liouville–von Neumann equation, such as the influences of radiofrequency (RF) pulses and free precession periods, are considered next.


## 11. Radiofrequency Pulses

A pulse sequence that can be used to reveal the split-component NMR spectra of quadrupolar nuclei in anisotropic surroundings is shown in Figure2. It consists of a typical pulse-acquire block, i.e., a 90x RF pulse of duration τ90 followed by signal acquisition (-y digitizer phase) during a time τacq. In the following sections, we will show how the residual quadrupolar coupling constant (RQC) invoked by the anisotropic environment of a stretched hydrogel is encoded into the state of the system during the signal acquisition time period and how the information subsequently becomes manifested in the resulting NMR spectrum.Figure 2
Pulse sequence used to obtain split-component NMR spectra of materials containing the quadrupolar nuclei2H, 23Na, and 133Cs as guests within anisotropic (stretched hydrogel) environments.The application of RF pulses shuffles the population distribution of nuclear spin states. A 90° RF pulse equalizes all state populations, while a 180° RF pulse inverts the spin state population hierarchy. The Hamiltonian for an RF pulseH^RF is as follows:(100)H^RF=ω1I^xcosω0t+ϕRF+I^ysinω0t+ϕRF,where ω1=γB1 is the amplitude of the RF field, expressed as a frequency, ω0=γB0 is the nuclear Larmor frequency, and ϕRF is the phase of the RF field relative to a specified reference. B1 corresponds to the magnitude of the magnetic field induced by a current flowing through a transverse coil used for signal detection. The complicated time dependence of equation (100) can be removed by transforming into a frame which precesses at a frequency ω0 for the spin(s) of interest, i.e., the RF pulse frequency for the spin(s), about the axis of the static magnetic field B0, which is assumed to be the z-axis of the laboratory frame [25]. In this rotating frame, the RF pulse Hamiltonian has the following form:(101)H^RF′=ω1I^xcosϕRF+I^ysinϕRF,where ′denotes the rotating frame RF pulse Hamiltonian. This frame transformation is achieved by assuming that the RF pulse is applied on resonance, i.e., the frequency of the RF pulse matches the nuclear Larmor frequency for the spin(s) of interest. Moreover, the oscillating magnetic field, by convention, is applied along the x-axis of a Cartesian frame of reference (ϕRF = 0), further simplifying the RF Hamiltonian to(102)H^RF′=ω1I^x.It should also be noted that the properties (shape over time, intensity, and duration) of the applied pulsed RF field are important for stimulating nuclear transitions, e.g., the modulation frequency of the applied oscillating field should match the nuclear Larmor frequency at the resonance condition for the spin(s) of interest.Consider the initial state of a spin system to be the magnetization present at thermal equilibriumρ^0=I^z. The identity operation has been dropped from the spin density operator for simplicity since all operators commute with the identity operator and commutation is linear. Since the operator I^x does not commute with the operator I^z, the cyclic commutation operations given in equation (26) apply. The evolution of the initial spin state under the effects of the rotating frame RF Hamiltonian is described by the Liouville–von Neumann equation:(103)ρ^1=−iℏH^^RF′ρ^0,=−iℏH^RF′,ρ^0,=−iℏH^RF′ρ^0−ρ^0H^RF′,where ρ^1 is the spin density operator after the application of the rotating frame RF Hamiltonian. And, we have introduced the nomenclatural convenience of the commutation superoperator (double caret) which performs the commutation operation, in this case, of H^RF′ with ρ^0. We dropped the time argument (t) to simplify the expression. To further simplify matters, evolution of the density operator ρ^0 under incoherent relaxation processes has also been neglected. The true form of ρ^1 is revealed by applying a rotation sandwich [25] to ρ^0, where the sandwich is constructed from rotation operators corresponding to the chosen axis of rotation (see equation (19)) which uses the property that rotation matrices can be easily obtained from the exponential of a Hamiltonian operator, and is(104)ρ^1=cosω1τβI^z−sinω1τβI^y,where τβ is the time over which the rotation of the magnetization vector by the RF pulse takes place through an angle β. The equation is simplified in the case of a β=90° RF pulse to(105)ρ^1=−I^yforτ90=2π4ω1.By using operator relationships that are analogous to those given in equations (43)–(45), we can express ρ^1 in terms of the raising and lowering operators (in analogy with equation (43)):(106)I^x=I^++I^−2,(107)I^y=I^+−I^−2i,thus enabling ρ^1 in equation (105) to be rewritten as(108)ρ^1=−12iI^+−I^−.By applying the RF field perpendicular to the static magnetic fieldB0, spin transitions with energy intensities that are proportional to the loss or gain of angular momentum are generated. Equation (108) suggests that by applying the pulsed magnetic field parallel to B0, no observable NMR transitions would arise since the operator I^z commutes with itself.


## 12. Free Evolution

An essential part of many NMR experiments is free evolution of spin states that are perturbed from equilibrium by RF pulses. Suppose that an ensemble of nuclear spins is described by a spin density operatorρ^1 at an arbitrary time t1. The state of the spin system could be the initial magnetization available at thermal equilibrium or a delicately prepared population distribution amongst all spin states established by using intricately designed pulse sequences. The spin system is allowed to evolve freely, without interference from RF pulses or changes in magnetic field, for a period τ to a time point t2. The density operator specifying the properties of the spin ensemble at the time t2 is denoted by ρ^2 and can be readily predicted by using the theory we have (and will develop further) here. The case of free evolution of a spin system at thermal equilibrium is rather uninteresting, so we will now consider the state of a spin system ρ^1 that is created by applying a 90° RF pulse.It is worth noting that in many experiments, the influence of free evolution under a chemical shift Hamiltonian may additionally be present, for example, if an experiment includes the use of an offset frequency. For simplicity, we have neglected this interaction and we also assume that the RF pulse is on resonance with the spin system of interest at all times.The free evolution of the density operator is described by using the operator formalism presented above (equation (103)). As in the case of applied RF fields, the relationship between the density operator at time points t1 and t2 is again governed by evolution that is described by the Liouville–von Neumann equation:(109)ρ^2=U^^21ρ^1,=U^21ρ^1U^21†.This remains as the solution of the equation of motion for the spin ensemble, but with the rotating frame RF Hamiltonian commutation superoperatorH^^RF′ replaced by the free evolution propagation superoperator U^^21τ. The operator U^21 is defined in equation (99), and the matrix representation of U^21τ takes the exponential form of the nuclear spin Hamiltonian. Therefore, it is evident that the Hamiltonian describes the free evolution of the spin system. The presence of a free evolution period is of paramount importance to the prediction of the relative intensities of multicomponent NMR spectra, since the effects of RQCs are encoded within the density operator.The form of the propagation operator allows the deterministic prediction of a target spin density operator at some time in the future, given an initial set of starting conditions. In the case of free evolution of a spin system under the influence of a quadrupolar HamiltonianH^Q, the resulting form of ρ^2 is(110)ρ^2=−12ie−iH^Qτ/ℏI^+−I^−eiH^Qτ/ℏ.At this point, it is now evident that there are terms in the density operator which are proportional to the quadrupolar coupling constant. However, we have not currently explored the full form of the quadrupolar Hamiltonian, and so as of yet, we do not know the commutation properties ofH^Q with the operators of complex magnetization. Furthermore, taking the effects of anisotropic interactions between nuclear spins and their surrounding environment into account is not straightforward. To understand how this key phenomenon works in the case of RQCs resulting from anisotropic surroundings and to decode the problem posed in the title of this article, we first need to understand the formalism of the quadrupolar Hamiltonian.


## 13. Quadrupolar NMR

Quadrupolar nuclei constitute ∼74% of the magnetic nuclei in the isotopes of all the elements of the periodic table. Therefore, comprehending the underlying nature of the quadrupolar interaction is of great importance in NMR spectroscopy. A key property of a quadrupolar nucleus is its electric quadrupole moment, which interacts with local electrical field gradients (EFGs) at the site of the nucleus. The EFGs can be represented by a tensor with three principal axis elements (see Figure3(a)). The sum of the principal axis elements is, by definition, zero: Vxx+Vyy+Vzz=0, i.e., it is a traceless tensor. The largest principal value of the EFG tensor Vzz is typically specified as eq which is the EFG at the quadrupolar nucleus [27]. The parameter ηQ describes the biaxiality of the quadrupolar interaction: ηQ=Vxx−Vyy/Vzz. In the following analysis, the quadrupolar interaction assumes a biaxiality of zero and hence is a uniaxial EFG tensor.Figure 3
(a) Ellipsoid representing a uniaxial quadrupole tensor withVxx≤Vyy≤Vzz. Schematic representations of quadrupolar nuclei with uniaxial electrical field gradient (EFG) tensors dissolved in (b) an isotropic hydrogel and (c) an anisotropic hydrogel. Stretching the hydrogel along a particular axis creates the anisotropic environment. The extent of alignment of the quadrupolar nuclei with the polymer materials of the hydrogel is dependent on the extent of uniform stretching.
(a)(b)(c)The rank-2 spherical tensor operators (STOs) are useful in constructing the quadrupolar Hamiltonian for a spin with quantum numberI > 1/2 (see Figure 4). The physical meaning of an STO with rank l and order m is based on its behaviour under rotation in three-dimensional (3D) space (see Appendix A [28]). For the case of an ensemble of noninteracting quadrupolar nuclei with spin I > 1/2, the nuclear spin Hamiltonian is written as(111)ℋ^Q=ωQ6T^20,where ωQ is the quadrupolar coupling constant (a real-valued constant) and T^20 is an STO characterizing the spin-related portion of the quadrupolar interaction. The quadrupolar Hamiltonian given in equation (111) is written to first-order and is equivalent to the secular approximation [25]. In this case, the quadrupolar Hamiltonian is only sufficient to describe the quadrupolar interaction when ωQ is small, i.e., orders of magnitude smaller than the nuclear Larmor frequency, which is always assumed in this article, and is very relevant for the case of RQCs.Figure 4
Visual representation of the spherical tensor operators (STOs) used to construct the nuclear spin Hamiltonian of a quadrupolar nucleus. The STOsT^lm are specified by their rank l and order m (see Notebook 3 for details regarding how to reproduce this figure).The quantitiesωQ and T^20 are defined as follows [25, 27]:(112)ωQ=3eQV¯zz2I2I−1ℏ,T^20=3I^z2−II+1I^6,where Q is the electric quadrupole moment of the nucleus [27], V¯zz is the zz-component of molecular length-scale EFG tensor averaged over internal motions and overall rotations, and I^ is the identity operator.In isotropic liquids, the quadrupolar coupling (and first-order quadrupolar Hamiltonian) averages to zero due to molecular tumbling in solution. The motionally averaged quadrupolar Hamiltonian has no influence on the spectral position of NMR peaks in isotropic environments in the case of a small residual quadrupolar coupling constant.In anisotropic media, such as stretched hydrogels, liquid crystals, and strained lattice structures [29], ωQ does not vanish under molecular reorientation and the alignment of the EFG tensor principal axes with respect to the static magnetic field B0 can often lead to a sophisticated relationship between molecular orientation and the resulting quadrupolar interaction. This complicated dependence can be manifest in the NMR spectrum of a sample of quadrupolar nuclei as a complicated lineshape of split-multiplet structure with intensities related to the nuclear spin quantum number I (see above). In general, the zz-component of the EFG tensor is given by a weighted average over all molecular orientations:(113)V¯zzϑ=∫VzzϑPϑdϑ,where Vzzϑ is the EFG for a molecular orientation ϑ and Pϑ is the probability that a molecular orientation is occupied. After creating an anisotropic phase, as in the case achieved by stretching a hydrogel, the molecular orientation dependent EFG V¯zzϑ becomes nonzero. Consequently, ωQ is also finite (an effective value which is not determined by molecular tumbling) and is referred to as the RQC constant throughout this article. We denote this residual quantity by ω¯Q in the remainder of the article. It is worth stressing at this point that it is the average value of the residual quadrupolar coupling constant ω¯Q that determines the size of the splitting between the multiplet components in the NMR spectra of quadrupolar nuclei in anisotropic environments. However, the peak linewidths are predominantly governed by relaxation processes, and in this case, it is the nonaveraged value of the residual quadrupolar coupling constant ωQ which is important. We leave further discussion of this point to another forum, and now focus on the main topic of this article.


## 14. Interpreting NMR Spectra

Having previously written the matrix representation ofI^− (since elsewhere we have asserted that this is the operator whose expectation value returns the spectrum we seek), all we need to do now is to calculate its expectation value.The frequency and amplitude of transitions between populated states are encoded within the spin density operator and are inferred by determining the amount of the observable operatorI^−that is present in the spin density operator ρ^, i.e., the projection of the lowering operator I^− onto the spin density operator ρ^. The reason for observing I^−, rather than I^x or I^y, allows for the detection of complex magnetization (see equations (105) and (106)), which encodes phase sensitive information.The projection simply selects the relevant matrix elements ofρ^ linking the basis states. Projections of matrix components with a common index (the diagonal elements) correspond to population operators, i.e., the probability of finding a spin in a pure state. The shift operators between orthogonal states denote coherences, if the expectation value is nonzero. This will describe nuclear transitions between the pure states and hence detectable magnetization that will be manipulated to produce the emergent NMR spectrum.With this in mind, the expectation value of the operatorI^− is determined by invoking equation (92):(114)I^−=Trρ^2I^−,where I^− is the observed outcome induced by the detection of the NMR signal. This result infers which elements of the spin density operator ρ^2 are observable. Off-diagonal elements of the density matrix contribute to the NMR signal, whereas diagonal elements of ρ^2 only contribute to longitudinal magnetization. The definition of ρ^2 in equation (110) can be used to expand equation (114), and as a result, the final outcome is(115)I^−∝TrI^+I^−−I^−I^−.The termI^+I^− produces diagonal matrix elements, with these components featuring in the NMR spectrum. This is the reason why the I^− operator is of great importance to solving a large number of questions in magnetic resonance spectroscopy. The term I^−I^− leads to off-diagonal matrix elements (the matrix representation of the I^−I^− term is traceless) that do not feature directly in the NMR spectrum, and as such simply implementing the trace (Tr) operation is sufficient to determine the features of the acquired NMR spectrum.The number of terms produced by equation (115), and the proportionality of the relationship between each term, must now be explained: this problem can be understood by considering the matrix dimensions of the lowering operator I^−. The number of terms is simply proportional to the spin angular momentum quantum number I, with the number of terms being 2I. These terms represent peaks in the acquired NMR spectrum, with each peak separated by the value of the RQC constant (see below). The amplitude coefficients of each term pose a significantly more difficult question. However, as hinted at previously, the anticipated proportionalities are controlled by the matrix structure of the lowering operator I^−.AMathematica calculation (see Notebook 2) shows that the relative coherence amplitudes of the different transitions between nuclear spin states are(116)2Ifor±I⟺±I∓1,(117)22I−1for±I∓1⟺±I∓2,(118)6I−1for±I∓2⟺±I∓3,(119)22I−3for±I∓3⟺±I∓4,where ⇔ represents a transition between two spin states separated by a single unit of spin angular momentum h. Of course, the transition amplitudes and their related matrix elements depend on the value of I since, e.g., a spin I = 1 nucleus will not have a transition between ±I∓3 and ±I∓4 states. Nevertheless, these transition amplitudes are identical to the matrix elements of the ladder operator with spin quantum number I (see equations (78) and (79)). It is at this point in our analysis that we have uncovered the root of the relative intensities of NMR spectral multiplets for quadrupolar nuclei in anisotropic media. The transition amplitudes given above are what we need for the remainder of this article. In principle, this sequence of values can continue indefinitely and different scaling factors can also be applied in successive calculations; however, the relative peak intensities do not waiver.In the context of this article, the focus is on the RQC remaining for a quadrupolar nucleus embedded within a partially aligned medium. Specifically, quadrupolar nuclei of salts in stretched hydrogels, which lead to the predictable intensities of split-multiplet NMR spectra (see Figure1). In the sections below, the mechanism of the degeneracy breaking interaction is included, and the manifestation of this phenomenon in the resulting NMR spectrum is scrutinized by using a product operator analysis for a number of different cases of the spin angular momentum quantum number I. Before showing specific examples, it is worth taking a moment to consider the properties of a stretched hydrogel.


## 15. Stretched Hydrogels

Hydrogel samples containing dissolved salts of quadrupolar nuclei, or of molecules with atoms with quadrupolar nuclei, usually display a single NMR peak, or resonance, that has a Lorentzian lineshape [20, 25]. Stretching the hydrogels breaks the degeneracy of the nuclear Zeeman transitions (see Figure 3), since the created anisotropic environment possesses molecular-length-scale EFGs [1, 17, 30–33]. The degeneracy breaking is known to originate from the quadrupolar nuclei possessing an electric quadrupole moment whereby interactions occur with the EFG tensor of the partially aligned polymer molecules of the hydrogel and is also related to a change in alignment tensors due to stretched hydrogel pore shapes. Prolate pore shapes lead to favoured molecular alignments and consequently preferential orientations of rotational diffusion [34]. Predictable multiplet NMR spectra with characteristic relative intensities of their components arise; the number and their relative intensities are specific to the quantum spin number I of the involved nucleus [33].


## 16. Previous Works

In a previous paper [35], an explanation for the relative intensities of the multiplet components was based on as little quantum mechanical theory as possible, and it relied on three basic ideas: (i) NMR-active atomic nuclei can occupy a finite number of discrete energy levels of which the total number is specific to the nuclide; (ii) the relative proportion or distribution of the nuclei in states of different energies is specified by a Boltzmann distribution [19]; and (iii) the NMR spectral outcome from the nuclear system, after equalizing the number of nuclei in each energy level, depends on the net number of swappings (“tumbling down the energy cascade”) that occurs between successive energy levels, on re-establishing the Boltzmann equilibrium state.While the explanation for the spectral features noted above is relevant in some pedagogical contexts, the exact theory of quantum mechanics is the final arbiter for a quantitative analysis because of the extent of its applicability across all forms of spectroscopy. Eliav [36] has recently provided an explanation of the phenomenon at hand based on fictitious spin-1/2 operators, which introduces a higher level mathematical construct that takes us away from our simpler pedagogical route and is not required here. However, in the present paper, we expand on the theory first given in [1, 17, 33, 35, 36] and posit that this theory constitutes a logically complete explanation for the relative intensities of the components for numerous examples of quadrupolar-split multiplets.


## 17. The Case ofI = 1

We have now assembled all the tools necessary to predict the relative intensities of split-component NMR spectra of solutes containing quadrupolar nuclei embedded within stretched hydrogels. The simplest case is the quadrupolar nucleus2H, which in isotropic solution displays a single resonance line in the 2H NMR spectrum. 2H nuclei in water have previously been employed to investigate the alignment of cartilage [37] and other oriented media [38].The elements in the matrix representation of the quadrupolar Hamiltonian depends on the chosen basis set and in this case can be expressed in the basis of kets of the angular momentum operator:1,+1, 1,0, 1,−1. Under these circumstances, the matrix representation of ℋ^Q for an I = 1 nucleus is(120)ℋ^Q=ω¯Q6000−ω¯Q3000ω¯Q6.The matrix representation of the Hamiltonian is diagonal since theT^20 operator is likewise diagonal:(121)T^20=16000−2300016.The initial state of a spin system for the majority of conventional NMR experiments is the usable magnetization present at thermal equilibrium. In such instances, the spin density operator is proportional to thez-angular momentum operator, namely, ρ^0=I^z. Since the longitudinal angular momentum operator corresponds to differences in populations between states, I^z is strictly diagonal and has the following matrix form:(122)I^z=10000000−1.As discussed in the section on RF pulses, the application of an RF field along thex-axis of a Cartesian frame of reference rotates the direction of the ensemble angular momentum towards the -y-axis. At this point the spin density operator is given by ρ^1=−I^y. The -y-angular momentum operator corresponds to coherences between all states [20, 25]. All spin state populations have also been equalized by the application of the RF field, and hence the diagonal components of −I^y are zero. The matrix representation of −I^y in the same basis is as follows:(123)−I^y=0i20−i20i20−i20.It will soon become apparent that the FID is a key part of the experiment. The spin system evolves during the signal acquisition timeτacq under the free evolution propagation operator U^21, which acts to transform ρ^1 into ρ^2. The rules as to how the propagation operator behaves are described in the section on free evolution. The propagation operator U^21 is defined as follows:(124)U^21=e−iℋ^Qτacq/ℏ.Equation (124) has the same form as equation (99). The propagation operator U^21 has the following matrix form:(125)U^21=e−iω¯Qτacq/6ℏ000eiω¯Qτacq/3ℏ000e−iω¯Qτacq/6ℏ,since U^21 is simply the matrix exponential of the quadrupolar Hamiltonian. ρ^1 is propagated forward in time by U^21 to a state ρ^2, which has the following form:(126)ρ^2=U^21ρ^1U^21†,=−e−iℋ^Qτacq/ℏI^yeiℋ^Qτacq/ℏ,where U^21† is the Hermitian adjoint (introduced in equation (109)) of the free evolution propagation operator. Equation (126) has the same form as equation (98). The matrix form of ρ^2 is written as follows:(127)ρ^2=0i2e−iω¯Qτacq/2ℏ0−i2eiω¯Qτacq/2ℏ0i2eiω¯Qτacq/2ℏ0−i2e−iω¯Qτacq/2ℏ0.The components ofρ^1 have acquired a frequency factor that is proportional to the size of the RQC. The free evolution period has successfully encoded information regarding the RQC into the spin density matrix ρ^2. At this point, the RQC can have a remarkable influence on the acquired NMR spectrum. The components of the resulting NMR spectrum are determined by equation (114), i.e., the trace (Tr) of the product between ρ^2 and the NMR observable operator I^−, which has the following matrix representation (see equation (116)):(128)I^−=0002I0002I0,=000200020.Detectable elements of the NMR signal, as inferred from equation (115), are selected by employing the trace operation:(129)Trρ^2I^−=2e−iω¯Qτacq/2ℏ+2e−iω¯Qτacq/2ℏ.The two terms in equation (129) encode two vital pieces of information with regard to the form of the acquired NMR spectrum: (i) the exponential factor determines the resonance position of the peak; and (ii) the coefficient gives the relative intensity or amplitude of the allowed transition. The two transitions would become degenerate in the absence of an RQC, and the overall intensity of the detected NMR peak would be 4. Such a scenario is present for the case of quadrupolar nuclei in a relaxed gel or in isotropic solvents. For an I = 1 quadrupolar nucleus, such as 2H, the resonance positions in the NMR spectrum appear at ±ω¯Q/2, i.e., the two peaks are separated by the size of the RQC constant. In this case, the intensities of the two peaks present in the 2H NMR spectrum are shared equally in a relative 2 : 2 (1 : 1) ratio. The relative intensity ratios are in agreement with equation (116). This result is also in agreement with previously predicted ratios for the relative signal intensities of 2H nuclei containing materials in stretched hydrogels [35]. A simulated spectrum of a 2H (I = 1) nucleus in an anisotropic environment is shown in Figure 5. The NMR spectrum was simulated by using the Mathematica-based NMR software package SpinDynamica [39] (see Appendix B for a snapshot of the SpinDynamica code used in this simulation). The simulation uses an RQC constant of ω¯Q/2π=360 Hz (see Notebook 4 for the simulation protocol employed to generate Figure 5 and the product operator analysis which yields equation (129)). Prediction of the relative split-component intensities for I = 1 nuclei is simple but we shall see that the picture becomes more complex when calculating the NMR spectra of higher spin nuclei. In the following sections, we address the experimentally relevant [3] cases of I = 3/2 and 7/2 nuclei in anisotropic media.Figure 5
Simulated NMR spectrum of a solute containing the quadrupolar nucleus2H (I = 1) dissolved in an anisotropic environment, such as a stretched hydrogel. The splitting between the peaks can be controlled and is a function of the anisotropy, i.e., the extent of stretching the hydrogel. In this case, the splitting between the peaks is 360 Hz, which corresponds to a RQC constant of ω¯Q/2π=360 Hz. The ratio of the relative peak intensities is: 2 : 2 (1 : 1). The NMR spectrum was simulated by using SpinDynamica [39]. Lorentzian line broadening (full-width at half-maximum height = 40 Hz) was used in the simulation.


## 18. The Case ofI = 3/2

For spinI = 3/2 nuclei, the additional NMR spectral complexity arises from the increased number of basis states. However, as shown later, this neatly leads to the consequent realization of an additional peak in the NMR spectrum and an alternative ratio of the relative peak intensities for the split-multiplet components. Moreover, the larger matrix size can still be handled easily by current mathematical software packages such as Mathematica.In this section, identical thought processes and simulation protocols are implemented to obtain the correct (experimentally observed) multiplet structure for a spinI = 3/2 nucleus hosted by an anisotropic medium. A commonly encountered I = 3/2 quadrupolar nucleus is 23Na, which is often used in biomedical experiments [40, 41]. The quadrupolar Hamiltonian matrix representation of a 23Na nucleus depends on the spin I = 3/2 basis:(130)ℋ^Q=ω¯Q20000−ω¯Q20000−ω¯Q20000ω¯Q2,where the basis kets are as follows: 3/2,+3/2, 3/2,+1/2, 3/2,−1/2, 3/2,−3/2. The strictly diagonal matrix representation of ℋ^Q is due to the diagonal nature of the T^20 operator in this basis. ℋ^Q is of a similar form to equation (120).Consider the magnetization available at thermal equilibrium. In this case, the starting spin density operator corresponds to longitudinal magnetizationρ^0=I^z, and has the following matrix form:(131)I^z=320000120000−120000−32.In words, the matrix representation ofI^z has finite elements along the diagonal, with all other elements being zeros. Consider the pulse sequence presented in Figure 2. An RF field applied with phase x in a Cartesian axis system will rotate the initial z-magnetization into the -y-axis. The matrix representation of ρ^1=−I^y (in the same basis) is as follows:(132)−I^y=0i32−i32000i00−i000i32−i320.The matrix representation of−I^y has zeros on the diagonal and nonzero entries either side of the diagonal. These matrix elements are associated with single quantum coherences [20, 25].The timeτacq allows the spin system to evolve and as a result encode information regarding the RQC into the spin density operator. This is the change of ρ^1 into ρ^2, which is achieved by the action of the free evolution propagation operator U^21 defined in equation (99). In this case, the propagation operator U^21 has the following matrix form:(133)U^21=e−iω¯Qτacq/2ℏ0000eiω¯Qτacq/2ℏ0000eiω¯Qτacq/2ℏ0000e−iω¯Qτacq/2ℏ,where ρ^1 is propagated forward in time (described by equation (98)) and leads to the result ρ^2, which has the following matrix representation:(134)ρ^2=0i32e−iω¯Qτacq/2ℏ00−i32eiω¯Qτacq/2ℏ0i00−i0i32eiω¯Qτacq/2ℏ00−i32e−iω¯Qτacq/2ℏ0.The matrix form of this outcome contains exponential terms, which are related to the evolution of the spin system in the presence of the RQC. The coherences expressed in equation (134) are converted to observable magnetization during signal acquisition. The observable state of the spin system is derived by calculating the expectation value of the I^− lowering operator from the spin density operator during signal acquisition, i.e., acting with a measurement upon the form of the spin system in the state ρ^2 (see equation (114)). The matrix representation of I^−, in the same basis, is given below (see equations (116) and (117)):(135)I^−=00002I000022I−100002I0,=0000300002000030.The trace of this product leads to the detectable elements of the NMR signal:(136)Trρ^2I^−=3eiω¯Qτacq/ℏ+4+3e−iω¯Qτacq/ℏ.In the case of a spin-3/2 nucleus, there are three peaks in its NMR spectrum with the exponential factors in equation (136) indicating to which peak each intensity is associated. A simulated spectrum of a 23Na nucleus in an anisotropic environment is shown in Figure 6 [39]. The simulation uses an RQC constant of ω¯Q/2π = 180 Hz; it depicts a triplet multiplet structure. For the I = 3/2 quadrupolar nucleus, the resonance positions in the NMR spectrum appear at 0,±ω¯Q, i.e., the peaks are separated by the size of the RQC constant, as in the case of a spin-1 quadrupolar nucleus. Experimentally, the splitting between the peaks can be readily controlled as a function of the extent of gel stretching (see Figure 7). The size of the RQC increases with an increasing stretching factor. In this case, the intensities of the three peaks present in the 23Na NMR spectrum are shared in a relative 3 : 4 : 3 ratio, as predicted by the above theory (see equations (116) and (117)). This result is in agreement with previously observed ratios for the relative signal intensities of 23Na containing materials in stretched hydrogels [2]. If the RQC constant were equal to zero, as in the case of isotropic liquids, the NMR spectrum of the 23Na nucleus would be a single Lorentzian lineshape.Figure 6
Simulated NMR spectrum of a solute containing the quadrupolar nucleus23Na (I = 3/2) dissolved in an anisotropic environment, such as a stretched hydrogel. The splitting between the peaks can be controlled and is a function of the anisotropy, i.e., the extent of stretching the hydrogel. In this case, the splitting between the peaks is 180 Hz, which corresponds to a residual quadrupolar coupling (RQC) constant of ω¯Q/2π  = 180 Hz. The ratio of the relative peak intensities is 3 : 4 : 3. The NMR spectrum was simulated by using SpinDynamica [39]. Lorentzian line broadening (full-width at half-maximum height = 40 Hz) was used in the simulation.Figure 7
Experimental residual quadrupolar coupling (RQC) constants for (a)2H and (b) 23Na nuclei containing salts dissolved in 14% w/v bovine gelatin acquired at 25°C and 9.4 T as a function of the hydrogel stretching factor.
(a)(b)


## 19. The Case ofI = 7/2

The relative multiplet intensities predicted above were formulated by manipulating the spin density operator such that aπ/2 pulse followed by a free evolution period would lead to the manifestation of an RQC in the detected NMR spectrum. The formalism described above should be able to replicate the information in an NMR spectrum for any quadrupolar nucleus in a uniformly anisotropic environment. One such, perhaps slightly more tricky example, is the NMR spectrum of an I = 7/2 species hosted in an anisotropic medium. 133Cs+ ions, for which I = 7/2, have been used in NMR studies of cation transport across the membranes of red blood cells [4, 42, 43].In one final example, we now consider the much more complicated case of133Cs+ ions embedded in a hydrogel that is stretched. In this case, the complete basis set of kets for the single-nucleus system is 7/2,+7/2, 7/2,+5/2, 7/2,+3/2, 7/2,+1/2, 7/2,−1/2, 7/2,−3/2, 7/2,−5/2, and 7/2,−7/2. The quadrupolar Hamiltonian for the system contains an STO denoted by T^20 that is related to the angular momentum operators. It is this operator that gives ℋ^Q its diagonal matrix representation:(137)ℋ^Q=7ω¯Q200000000ω¯Q200000000−3ω¯Q200000000−5ω¯Q200000000−5ω¯Q200000000−3ω¯Q200000000ω¯Q2000000007ω¯Q2.The pulse sequence presented in Figure2 acts on the initial spin density operator ρ^0, which corresponds to thermal equilibrium magnetization,which is proportional to the z-angular momentum operator I^z, and has the following matrix form:(138)I^z=7200000000520000000032000000001200000000−1200000000−3200000000−5200000000−72.The evolution of the initial stateρ^0 under the effect of a rotating frame RF Hamiltonian, i.e., a π/2 pulse applied with x-phase, is described by the Liouville–von Neumann equation (see equation (103)). Solving the commutator equation leads to the matrix representation of ρ^1, which in this case is proportional to –I^y:(139)−I^y=0i72000000−i720i3000000−i30i152000000−i15202i000000−2i0i152000000−i1520i3000000−i30i72000000−i720.Equation (139) is of a similar form to equations (123) and (132).As discussed above, the free evolution delay is critical for realizing the NMR multiplet structure of quadrupolar nuclei in anisotropic surroundings. The evolution of the spin density operatorρ^1 is controlled by the propagation operator U^21, which is a function of the acquisition time τacq. The propagation operator U^21 is defined in equation (99); it has the following matrix form (in the same basis):(140)U^21=e−7iω¯Qτacq/2ℏ00000000e−iω¯Qτacq/2ℏ00000000e3iω¯Qτacq/2ℏ00000000e5iω¯Qτacq/2ℏ00000000e5iω¯Qτacq/2ℏ00000000e3iω¯Qτacq/2ℏ00000000e−iω¯Qτacq/2ℏ00000000e−7iω¯Qτacq/2ℏ,where ρ^2 is propagated forward in time by the free evolution operator U^21 to ρ^2, described by equation (98). ρ^2 provides the outcome for the state of the system after evolution under the quadrupolar Hamiltonian and has the following matrix form after τacq:(141)ρ^2=0i72e−3iω¯Qτacq/ℏ000000−i72e3iω¯Qτacq/ℏ0i3e−2iω¯Qτacq/ℏ000000−i3e2iω¯Qτacq/ℏ0i152e−iω¯Qτacq/ℏ000000−i152eiω¯Qτacq/ℏ02i000000−2i0i152eiω¯Qτacq/ℏ000000−i152e−iω¯Qτacq/ℏ0i3e2iω¯Qτacq/2ℏ000000−i3e−2iω¯Qτacq/2ℏ0i72e3iω¯Qτacq/2ℏ000000−i72e−3iω¯Qτacq/2ℏ0.The matrix elements of this spin density operatorρ^2 have acquired additional information during the free evolution interval, including an exponential term that is proportional to the magnitude of the RQC.Consider, at present, the successfully encoded information of the RQC as being contained within the spin density operatorρ^2. The multiplet structure of the desired NMR spectrum is unveiled by employing the NMR observable operator I^−, which has the following matrix representation:(142)I^−=000000002I0000000022I−1000000006I−10000000022I−3000000006I−10000000022I−1000000002I0,=00000000700000000230000000015000000004000000001500000000230000000070.The physical interpretation of a product betweenρ^2 and I^− is realized by simply asking the following question: How much of the lowering operator I^− is contained within the spin density operator ρ^2? The expectation value, i.e., the trace of the above mentioned product, described by equation (114), reveals the presence of the RQC within the acquired NMR spectrum, with the detectable components of the NMR spectrum, i.e., relative peak intensities and positions (offset frequencies), additionally extracted:(143)Trρ^2I^−=7e3iω¯Qτacq/ℏ+12e2iω¯Qτacq/ℏ+15eiω¯Qτacq/ℏ+16+15e−iω¯Qτacq/ℏ+12e−2iω¯Qτacq/ℏ+7e−3iω¯Qτacq/ℏ.In the case of a spin-7/2 nucleus, the resulting NMR spectrum has an iconic septet structure with a set of characterisitc relative intensities, as opposed to the case of isotropic solvents, which in most cases is a single Lorentzian resonance line. Figure8 shows a simulated NMR spectrum of 133Cs nuclei (I = 7/2) in an anisotropic environment [39]. The simulation used an RQC constant of ω¯Q/2π = 60 Hz. Such a spectrum arises from the presence of the seven terms in equation (143). The resonance positions of the peaks in the NMR spectrum of an I = 7/2 quadrupolar nucleus, such as 133Cs, constrained within an anisotropic medium, appear at 0,±ω¯Q,±2ω¯Q,±3ω¯Q. All peaks are separated by the magnitude of the RQC constant that is consistent with the spectra simulated for the cases of spin I = 1 and I = 3/2 nuclei in the above sections. The relative spectral intensities of the seven-peak multiplet NMR spectrum are given by the following enchanting seven term ratio: 7 : 12 : 15 : 16 : 15 : 12 : 7, in agreement with equations (116)–(119). This result is also in agreement with previously observed ratios for the relative signal intensities of 133Cs salts in stretched hydrogels [1] (see Figure 1). It is also worth pointing out that the NMR spectra in Figure 1 have previously been reproduced (dashed lines) using various simulation techniques.Figure 8
Simulated NMR spectrum of a solute containing the quadrupolar nucleus133Cs (I = 7/2) dissolved in an anisotropic environment, such as a stretched hydrogel. The splitting between the peaks can be controlled and is a function of the anisotropy, i.e., the extent of stretching the hydrogel. In this case, the splitting between the peaks is 60 Hz, which corresponds to a residual quadrupolar coupling (RQC) constant of ω¯Q/2π = 60 Hz. The ratio of the relative peak intensities is 7 : 12 : 15 : 16 : 15 : 12 : 7. The NMR spectrum was simulated by using SpinDynamica [39]. Lorentzian line broadening (full-width at half-maximum = 40 Hz) was used in the simulation.In this section, we demonstrated the power of quantum mechanical theory to predict the multiplet structure of NMR spectra from quadrupolar nuclei in the presence of an RQC, even in the complicated case of a nucleus with spinI = 7/2 [44]. Thus, our primary aim for this article has been realized.


## 20. Conclusions

A quantum mechanical analysis of quadrupolar-split NMR spectra has been presented. Basic mathematical tools were developed as part of guiding you, the reader, to the final result. The relative intensities and resonance positions (frequency offsets) of multicomponent NMR spectra were predicted for the cases of2H (I = 1), 23Na (I = 3/2), and 133Cs (I = 7/2) nuclei in uniformly anisotropic media. The results are in agreement with experimental NMR spectra of quadrupolar nuclei of inorganic salts in stretched hydrogels. It is worth noting that a similar reasoning would be adequate to describe the relative intensities of split-component NMR spectra for scalar coupled and residual dipolar coupled (RDC) systems, or “unfocussed” interactions. In the current analysis, the influence of RDCs was ignored, yet such interactions are regularly encountered in stretched hydrogel systems [3, 5–10].In Figure1, it is noted that there is a progressive broadening from the central transition outwards, with the outer components being approximately twice as broad as the central component. Although the components follow the predicted ratios for the relative integrals of the individual components, the qualitative difference in appearance between the NMR spectra presented in Figures 1 and 8 is due to the presence of a quadrupolar relaxation mechanism operating with the I = 7/2 nuclei. An extension to the quantum mechanical theory presented here would therefore be the inclusion of incoherent quadrupolar relaxation processes that are described by a Liouvillian superoperator approach.



---

*Source: 8890478-2021-02-19.xml*