An Introduction to the Holstein-Primakoff Transformation, with Applications in Magnetic Resonance
aa r X i v : . [ c ond - m a t . o t h e r] J un CONTENTS
I. Introduction 3II. Notations and definitions 7III. The Holstein-Primakoff Transformation 9A. The quantum harmonic oscillator 9B. The Holstein-Primakoff Transformation 16IV. Multispin magnetic resonance with the Holstein-Primakoff Transformation 23A. Index compression maps 23B. Eigendecomposition of Isotropic Multispin Hamiltonians 271. The isotropic multispin Hamiltonian in the absence of an external magneticfield 282. The isotropic multispin Hamiltonian in the presence of an external staticmagnetic field 513. Eigendecomposition of Liouvillians 57V. On Schwinger bosons 59VI. Final Remarks 63Acknowledgments 67References 671 n Introduction to the Holstein-Primakoff Transformation,with Applications in Magnetic Resonance
J.A. Gyamfi a Scuola Normale Superiore di Pisa, Piazza dei Cavalieri 7, 56126 Pisa, Italy.
Abstract
We have witnessed an impressive advancement in computer performance in the last couple ofdecades. One would therefore expect a trickling down of the benefits of this technological advance-ment to the borough of computational simulation of multispin magnetic resonance spectra, but thathas not been quite the case. Though some significant progress has been made, chiefly by Kuprovand collaborators, one cannot help but observe that there is still much to be done. In our view,the difficulties are not to be entirely ascribed to technology, but, rather, may mostly stem from theinadequacy of the conventional theoretical tools commonly used. We introduce in this paper a setof theoretical tools which can be employed in the description and efficient simulation of multispinmagnetic resonance spectra. The so-called Holstein-Primakoff transformation lies at the heart ofthese, and provides a very close connection to discrete mathematics (from graph theory to numbertheory).The aim of this paper is to provide a reasonably comprehensive and easy-to-understand intro-duction to the Holstein-Primakoff (HP) transformation (and related bosons) to researchers andstudents working in the field of magnetic resonance. We also focus on how through the use ofthe HP transformation, we can reformulate many challenging computing problems encountered inmultispin systems as enumerative combinatoric problems. This, one could say, is the HP transfor-mation’s primary forte. As a matter of illustration, our main concern here will be on the use of theHP bosons to characterize and eigendecompose a class of multispin Hamiltonians often employedin high-resolution magnetic resonance. a Email: jerryman.gyamfi@sns.it . INTRODUCTION Consider the deuterated hydroxymethyl radical ( • CH OD). It is well understood that theexperimental determination of the electron spin resonance (ESR) spectrum for the radicalinvolves the application of a static magnetic field B o and a microwave frequency field B ( t ) (or a succession of such pulses), perpendicular to B o . How these two magnetic fields areeffectively applied defines an experimental method. We could, for example, hold B o constantand vary the frequency of the irradiation field, or vice versa. In all cases, though, typicalexperimental procedures require subjecting the sample first to the static magnetic field B o , followed then – after an interval of time ∆ t sufficient enough to allow the system toconsiderably approach its equilibrium state – by the irradiation field B ( t ) .Suppose we want to simulate the ESR spectrum for the • CH OD radical. Like all quantummechanical calculations, we start with a suitable Hamiltonian. In this case, the Hamiltonianwill consist of that for the isolated radical, plus its interaction with the magnetic fields. Thefull Hamiltonian for the isolated radical alone (which takes into account the core electrons)is way too complicated to deal with. But, fortunately, experience has taught us that the spinHamiltonian approach, though phenomenological by nature, suffices for the task. With thisapproach, we ignore the core electrons altogether and take into account only the interactionbetween: 1) the unpaired electron and the external magnetic fields, 2) the magnetic nuclei(i.e. nuclei with nonzero spin quantum numbers) and the external magnetic fields and 3)the unpaired electron and the magnetic nuclei. Hence, in the case of • CH OD we only needto consider four particles with nonzero spin quantum numbers: the unpaired electron (spin- / ), the two methyl hydrogen nuclei (each of spin- / ) and the deuterium nucleus (spin- )– where we have obviously assumed that the carbon and oxygen nuclei are of the isotopes C and O, respectively, which are both zero spin nuclei. Let us assume that all couplingsare scalar, and that the spin-spin couplings are simply of the type T ii ′ ˆ J ˆ J ˆ J i · ˆ J ˆ J ˆ J i ′ [1], where T ii ′ isthe coupling constant and ˆ J ˆ J ˆ J i (= ˆ J xi e x + ˆ J yi e y + ˆ J zi e z ) is the spin vector operator for spin i .In the absence of the external magnetic fields, the spin Hamiltonian for the isolated radical, ˆ H spin − spin , is given by: ˆ H spin − spin = X i>i ′ T i,i ′ ˆ J ˆ J ˆ J i · ˆ J ˆ J ˆ J i ′ = X i>i ′ T i,i ′ (cid:16) ˆ J xi ˆ J xi ′ + ˆ J yi ˆ J yi ′ + ˆ J zi ˆ J zi ′ (cid:17) . (1)In the presence of the external static magnetic field B o , we need to add to the Hamiltonian3 H spin − spin an extra term which accounts for the interaction between the spins and B o . TheHamiltonian for the spins then becomes ˆ H o , where ˆ H o := − X i γ i ˆ J ˆ J ˆ J i · B o + ˆ H spin − spin (2)where γ i is the gyromagnetic ratio for the i − th spin. (We adopt in this article the conventionaccording to which the gyromagnetic ratio carries a sign; for example, it is negative for theelectron, and positive for the proton and most atomic nuclei). If we let the system approachits equilibrium state near enough, under the evolution of the Hamiltonian ˆ H o , then thecorresponding density matrix operator, ˆ ρ o , of the spin system prior to the application of theirradiation field B ( t ) could be taken – without loss of generality – to be of the form: ˆ ρ o = e − β ˆ H o Z (3)where β = k B T ( k B being the Boltzmann constant and T the absolute temperature), and Z is the quantum statistical partition function defined as: Z := Tr h e − β ˆ H o i . (4)The calculation of the expectation value of observables (thus, the magnetic resonance spec-trum by extrapolation) related to the spin system after the application of B ( t ) is highlydependent on ˆ ρ o . But given that the operator ˆ ρ o commutes with ˆ H o – the important im-plication being that they both share the same set of eigenvectors – the eigenvectors andeigenvalues of ˆ H o thus assume a fundamental role in the computational endeavor to sim-ulate the magnetic resonance spectrum of a spin system – because knowing them greatlysimplifies the computation.One simple way to determine the eigendecomposition of the operator ˆ H o is to take itsmatrix representation in some spin basis (the uncoupled spin representation being the mostcommon one) and diagonalize the whole matrix tout court (let us call this the "tout courtapproach"). But there is a catch here: as it is well known, the dimension of the Hilbertspace for a system of nonzero spins grows exponentially with the number of these spins.Thus, the dimension of the matrix representation of ˆ H o also inflates exponentially withthe number of nonzero spins. (For example, a system which consists of N spins- / hasa spin Hilbert space of dimension N .) The practical consequence of this – as far as theart of simulation is concerned – is that the computational cost of the eigendecomposition4lso scales exponentially with the number of spins, which in turn translates into far longercomputing time. There is also the problem of computing memory to mention: the finitememory available for use on computing devices means that the dimension of ˆ H o cannotexceed a given threshold inherent to the device. All these fall under what has been dubbedas the "curse of dimensionality". The mentioned inconveniences, together with the numericalinaccuracies one easily encounters when numerically diagonalizing very large matrices, castthe tout court approach as practically unsuitable for normal systems with more than aboutfifteen nonzero spins.Symmetry and arguments from the theory of angular momentum can also be used toreduce the computational complexity of eigendecomposing spin Hamiltonians of the form of ˆ H o . These arguments certainly depend on the symmetry properties of the spin Hamiltonianand may not be enough to alleviate the curse of dimensionality as we shall see later.In this article we show how one can further reduce the eigendecomposition computationalcost using the so-called Holstein-Primakoff (HP) representation, first introduced by TheodoreHolstein and Henry Primakoff in 1940 to study spin waves in ferromagnets[2]. Simply put,the HP representation is just another way of representing spin states and operators. So far,its application has been mainly confined to the spin wave theory in solid state physics but itis becoming increasingly clear that it has far reaching applications in the theory of angularmomentum and group theory than previously thought[3]. The fundamental role played byangular momentum theory in magnetic resonance theory cannot be overemphasized; thus,it is not surprising that the HP representation also finds important applications in thelatter. Our scope here is to illustrate some of these applications, namely, in relation to theeigendecomposition of spin Hamiltonians of the type ˆ H o , Eq. (2). If we take the • CH ODradical for example, the spin Hilbert space is of dimension 24. The tout court approach willhave us diagonalize a matrix of dimension 24 to find the eigenvectors and eigenvalues of theradical, but in the HP representation the problem reduces to effectively diagonalizing justtwo matrices: one of dimension 4 and the other of dimension 7. As we shall see, besidesenabling us to predict how many smaller matrices we will have to diagonalize and theirrespective dimensions, the HP transformation approach can also be used to infer the beststrategy to construct these matrices in a computationally efficient manner.But more importantly, our scope here is to introduce this type of spin representationto researchers and students – be they experimentalists or theorists – working in the field5f magnetic resonance. We therefore give here a very comprehensive introduction to thesubject before applying it. No advanced knowledge of quantum mechanics, group theory oreven some knowledge of number theory is required of the Reader. The scope of the paperis therefore twofold: 1) serve as a hard-to-come-by simple and comprehensive introductionto the HP transformation, and 2) show some novel applications of the latter in magneticresonance, with particular reference to multispin systems.What we have discussed above is the diagonalization approach to the simulation of mag-netic resonance spectra. An alternative is the propagation approach, where one propagatesthe density matrix ˆ ρ o , usually in Liouville space. One notable computational frameworkbased on this approach is the ingenious state space restriction (SSR) method proposed byKuprov and collaborators[4, 5]. The density matrix propagation approach is more appro-priate than the diagonalization approach when dealing with a large number of spins butwe limit ourselves in this article to the latter because it makes the exposition of the HPtransformation and the related mathematical apparatuses to be discussed more intelligible.Nevertheless, in the framework of HP transformation, the diagonalization approach in thecase of isotropic Hamiltonians can be pursued at a computational cost far less than the toutcourt approach, even though the scaling may still remain exponential with respect to thenumber of spins in some cases. The propagation approaches can all be reformulated in theframework of HP transformation, allowing the simulation of the resonance spectra at aneven lower computational cost.The paper is organized as follows: In §II we clarify the notations to be employed in thepaper, together with some definitions. The Holstein-Primakoff transformation is introducedin §III after we have revisited the quantum harmonic oscillator. This is to highlight some ofthe important features between the two. Then follows §IV, where we use the HP transforma-tion and other theoretical tools to be introduced there to, first, study the eigendecompositionproblem in relation to the isotropic spin Hamiltonian for an isolated multispin system; thatis then followed by a similar study on the same system but, this time, in the presence ofan external static magnetic field. In §V we introduce the Schwinger bosons, which is alsoanother type of spin representation but closely related to the HP transformation.6 I. NOTATIONS AND DEFINITIONS
To make the concepts we will discuss later on more intelligible, it is important we get ourterminologies and notations in order now. For example, we shall often speak of multisets [6].Unlike ordinary sets where each distinct element must appear only once, in a multisetdistinct elements can repeat any number of times. A set therefore can only tell us howmany distinct elements we have and their identities. For instance, say P the set of the firstfour prime numbers; then, P = { , , , } . Consider all the possible three digit (positive)integers we can create from the elements of P . The number "537" for example is a validone. By representing the digits composing these numbers as elements of a collection, wemay write "537" as " { , , } ", using the set notation. In the same spirit, { , , } , { , , } and { , , } are all valid. While { , , } is a set, { , , } , { , , } and { , , } are not,but are multisets. According to the criterion chosen to represent all possible three digitnumbers from the set P , it is clear that the multiset { , , } is different from { , , } .Order is therefore important here. Like a set, a multiset can also be ordered. If we areonly interested in the number of times a number repeats as a digit, then order is no longerimportant and so we may indicate "377" and "737" as { , , } or { , , } or { , , } , or –even more succinctly, using a customary notation in multiset theory – { , } , the exponent(or multiplicity ) here indicates the number of times an element repeats itself.Consider now a multiset A = { j , j , . . . , j i , . . . , j N } of N spins, where j i is the spinquantum number of the i − th spin. In the following, we shall reserve the Latin letter i to lower index A ’s elements when the latter is expressed in this laborious manner. A moreefficient way to indicate the same multiset (if we are not interested in order as we are now) isto specify all the distinct elements and their multiplicities. We shall use the Greek letter α toindex distinct elements of the multiset. Thus, if the spin multiset A has σ distinct elements,then we can express A also as A = { j N α α } , α = 1 , , . . . , σ , where N α is the multiplicity ofthe α − th distinct spin quantum number. Naturally, N = P σα =1 N α . It is crucial to pointout that two elements j i and j i ′ of A are distinct only if they correspond to different spinquantum numbers – other characteristics of the spins (like charge, magnetic moment etc.)are not of any merit whatsoever here. For example, the multiset A = {
12 9 , } indicates any aggregate of nine spin- / and three spin- ; the actual composition could consist of,for instance, 1) nine electrons and three protons, or 2) five N plus four muons – both7pin- / s – and three N nuclei (which are spin- )[7], etc.A very important distinction is due here. A system of N spins whereby j = j = . . . = j N = j is termed as univariate spin system (USS)[3]. The multiset representation of thesystem is then A = { j N } . In a USS, the N spins could be mutually different in regardsto mass, charge, magnetic moment, etc. On the other hand, a system of identical spins (IS)[3] is a univariate spin system whereby all the N spins share exactly the same intrinsicfundamental properties like mass and charge, and are indistinguishable from one anotherwhen placed under the same external conditions. Thus, an identical spin system is alsonecessarily USS, but a USS is in general not an IS.In addition to the concept of univariate and identical spin systems, we also have equivalent spins (ES). Consider a given multiset A = { j , j , . . . , j N } of spins. Say A ′ a submultiset of A , i.e. every element of A ′ is also an element of A . A ′ is said to be a system of equivalentspins if every element of A ′ couples to all other elements of A and any external field inthe same manner. Specifically, the elements of A ′ are said to be equivalent if one cannotdistinguish between them on the basis of their coupling tensors with other spins and externalfields. It mostly happens that equivalent spins are also identical, but in principle they donot need to be. The concept of equivalent nuclei in NMR, for example, is just a limit caseof equivalent spins. If we take the methyl radical • CH for example (assuming all threehydrogen nuclei are H and the carbon atom is C), the spin system A = { j , j , j , j } = {
12 4 } is clearly univariate; the three hydrogen nuclei are identical spins. The same trioof spins also constitute a collection of equivalent spins when the system’s Hamiltonian isinvariant under the operation of the point group C .Operators will be indicated with their usual hats while their matrix representations willbear none: for example, the matrix representation of the operator ˆ A will be simply indicatedas A .If j i is the spin quantum number of the i − th spin, therefore a scalar ( j i = , , , , , . . . ),then its corresponding spin vector operator will be indicated as ˆ J ˆ J ˆ J i , which, as we saw above,corresponds to the sum: ˆ J ˆ J ˆ J i = ˆ J xi e x + ˆ J yi e y + ˆ J zi e z , where ˆ J xi , ˆ J yi , ˆ J zi are the spin angularoperators along the respective axis and defined on the spin Hilbert space of spin i . In thisarticle, we will not try to indicate spin vector operators of electrons and nuclei with differentsymbols: they will all be indicated simply as ˆ J ˆ J ˆ J i , and the index i will serve the purpose ofrecording the specific identity of the spins based on our chosen choice of integer-labelling of8he latter. III. THE HOLSTEIN-PRIMAKOFF TRANSFORMATION
Before going on to discuss the Holstein-Primakoff transformation, it is useful we brieflyrecall here some salient features of the quantum harmonic oscillator. The reason why we doso is that the Holstein-Primakoff transformation, like many other techniques developed inmany-body quantum theory, can be traced back to the quantum harmonic oscillator. Theshort review we give below on the latter looks at the subject from a slightly different pointof view with respect to the traditional treatments, and focuses on those important featuresit shares with the HP transformation.
A. The quantum harmonic oscillator
The harmonic oscillator is perhaps the most important model in quantum mechanics.We are not going to belabor it here as it is extensively treated in almost every textbook onquantum mechanics. We only recall here that the Hamiltonian of a body of mass m tied toa spring able to move only in one direction (let us call it the x axis) is[8]: ˆ H = ˆ P x m + 12 k ˆ X (5)where ˆ P x is the linear momentum operator along the x − axis, ˆ X is the spatial displacementoperator from the equilibrium position, and k is the force constant.Since we shall be working in the position space, it is important to recall that the positioneigenkets {| x i} are such that: ˆ X | x i = x | x i (6)and they obey the orthogonality condition: h x ′ | x i = δ ( x ′ − x ) (7)where δ ( • ) is the Dirac delta function. If we multiply Eq. (7) by | x ′ i and integrate over theentire range of x ′ we get: Z dx ′ | x ′ i h x ′ | x i = Z dx ′ | x ′ i δ ( x ′ − x ) = | x i (8)9rom which we deduce that: Z dx ′ | x ′ i h x ′ | = ˆ I (9)where ˆ I is the identity operator. Eq. (9) is the completeness relation for the positioneigenkets.To study the behavior of the body quantum mechanically, what we can do is to solve thetime-independent Schrödinger’s equation for the body’s energy eigenkets, {| E i} : ˆ H | E i = E | E i ˆ P x m + 12 k ˆ X ! | E i = E | E i (10)where E is the energy of the body. Note that we may expand | E i in terms of the positioneigenkets: | E i = Z dx | x i h x | E i (11)where ( kh x | E ik dx ) can be interpreted as the probability to find the body in any of thepoints between x and x + dx given that its energy is E . The coefficient h x | E i =: ψ E ( x ) istherefore the probability amplitude, or wavefunction, in position space.If we multiply both sides of Eq. (10) from the left by the identity operator R dx | x i h x | ,we obtain: Z dx | x i (cid:18) m D x (cid:12)(cid:12)(cid:12) ˆ P x (cid:12)(cid:12)(cid:12) E E + 12 kx h x | E i (cid:19) = Z dx | x i E h x | E i . (12)Making use of the fact that D x (cid:12)(cid:12)(cid:12) ˆ P x (cid:12)(cid:12)(cid:12) E E = − ~ d dx h x | E i [9], Eq. (12) becomes: Z dx | x i (cid:18) − ~ m d dx + 12 kx (cid:19) ψ E ( x ) = Z dx | x i E ψ E ( x ) (13)that is, (cid:18) − ~ m d dx + 12 kx (cid:19) ψ E ( x ) = E ψ E ( x ) . (14)The solution to this eigenspectrum problem is well known. The energy eigenvalues are foundto be parameterized by n [8, 10]: E n = (cid:18) n + 12 (cid:19) ~ ω , n = 0 , , , , . . . (15)Thus, ψ E ( x ) = ψ ( n + ) ~ ω ( x ) , which may be simply written as ψ n ( x ) . Since there is a one-to-one correspondence between E and n (for fixed mass and force constant), we may henceforth10imply indicate | E i as | n i . Moreover, it is found that[8, 10]: ψ n ( x ) = 1 √ n n ! (cid:16) mωπ ~ (cid:17) / H n ( ζ ) e − ζ / (16)where ω ≡ p k/m , ζ ≡ x p mω/ ~ and H n ( z ) is the Hermite polynomial of order n in z .At this point we may rewrite Eq. (11) as: | n i = Z dx | x i ψ n ( x ) (17)where we have effected the substitutions | E i → | n i and h x | E i → ψ n ( x ) . It is evident fromEq. (17) that the wavefunction ψ n ( x ) is just the expansion coefficient when we expand | n i in terms of position eigenstates. Naturally, Eq. (17) is not the only possible expansion wecould think of. We could have equally expanded | n i in terms of the momentum eigenstates(along the x − axis), i.e. {| p x i} , where: ˆ P x | p x i = p x | p x i (18)and, h p ′ x | p x i = δ ( p ′ x − p x ) (19)from which follows the completeness relation: Z dp ′ x | p ′ x i h p ′ x | = ˆ I . (20)In fact, introducing the momentum space completeness relation into Eq. (17) we get: | n i = Z dp x | p x i (cid:18)Z dx h p x | x i ψ n ( x ) (cid:19) = Z dp x | p x i φ n ( p x ) (21)where φ n ( p x ) is the wavefunction in momentum space. In complete analogy to ψ n ( x ) , ( k φ n ( p x ) k dp x ) gives the probability of measuring the body’s momentum to be between p x and p x + dp x if its energy is fixed at E n . We briefly mention that because h p x | x i = √ π ~ exp (cid:0) − ip x x ~ (cid:1) [9], it follows that φ n ( p x ) = 1 √ π ~ Z dx exp (cid:18) − ip x x ~ (cid:19) ψ n ( x ) (22)which means φ n ( p x ) and ψ n ( x ) are related through the Fourier transform, as one wouldexpect.The gist of the above arguments is that depending on the basis in which we expand | n i ,we get different wavefunctions. The appropriate wavefunction to talk about depends on how11e intend to probe the quantum system: for example, if we intend to measure position givena fixed energy of the system, then { ψ n ( x ) } are the wavefunctions to go with; if, instead, wewant to measure momentum then we need the { φ n ( p x ) } . But in all these discussions, thenature of the ket | n i remains the same independent of whether we work in momentum orposition space. As a matter of course, one might then ask why don’t we just deal solelywith the energy eigenkets | n i without resorting to wavefunctions – and that is preciselythe quintessence of the so-called second quantization scheme (also known as the occupationnumber representation ). When we expand | n i in some basis and then solve for the relatedcoefficients, or wavefunctions, as we did above, Eq. (14), we speak then of first quantization .Thus, in first quantization, it is important to specify in which basis the expansion of the ket | n i is being done, normally in position or momentum space (but, in general, in the space ofeither one of a canonically conjugated pair of operators).How then do we go about doing quantum mechanics without explicitly talking of wave-functions? In other words, how does second quantization work? To see how it works, let usgo back to Eq. (10), which at this point, we may conveniently write as: ˆ H | n i = E n | n i . (23)Instead of trying to find an expression for | n i , we take it as it is – viz. {| n i} are eigenstatesof ˆ H , period. What we rather do is to rewrite the Hamiltonian ˆ H in terms of a set ofmutually commuting operators which have {| n i} as their eigenstates. It is worth noting thatthe eigenenergy E n , Eq. (15), depends only on n and not on neither the position ( x ) northe momentum ( p x ) of the body. Indeed, the mathematical expression for E n in terms of n gives important clues on how we may rewrite ˆ H in such a way that the new operators inthe latter do really have {| n i} as their eigenkets. If we combine Eqs. (23) and (15), we findthat, ˆ H | n i = (cid:18) n + 12 (cid:19) ~ ω | n i (24)which means that, ˆ H = (cid:18) ˆ A + 12 (cid:19) ~ ω (25)where ˆ A is, for now, an unknown operator with the following property: ˆ A | n i = n | n i . (26)12omparing Eq. (10) with Eq. (25), we find that: ˆ A = ωm ~ ˆ X + ˆ P ω m − ~ ωm ! . (27)Traditionally, this form of ˆ A is less preferred because it obscures some interesting insights.Rather, it is the factorized form: ˆ A = ˆ a † ˆ a (28)where[8], ˆ a † = r mω ~ (cid:18) ˆ X − imω ˆ P x (cid:19) (29a) ˆ a = r mω ~ (cid:18) ˆ X + imω ˆ P x (cid:19) (29b)that we prefer and use. The reason is that, despite the fact that the operators ˆ a † and ˆ a are not Hermitian, they allow for transition between the eigenstates – so they are extremelyuseful when discussing emission and absorption processes. As a matter of fact, one canderive that: ˆ a | n i = √ n | n − i , ˆ a † | n i = √ n + 1 | n + 1 i . (30)It is truly remarkable that the form ˆ H assumes in light of the operators ˆ a and ˆ a † , namely, ˆ H = (cid:18) ˆ a † ˆ a + 12 (cid:19) ~ ω (31)appears in many problems in quantum mechanics. Perhaps the most important is thequantization of the electromagnetic field, which was first done in the early years of quantumtheory. There, one sees that the Hamiltonian of the field becomes a sum over an infinitenumber of harmonic oscillators, each representing a specific mode of the field. And theoperators ˆ a † and ˆ a are tasked with increasing and decreasing by a quanta of energy theirrespective modes. Einstein had already discovered the photoelectric effect and the notion ofthe electromagnetic field as being composed of particles called photons was well establishedby then, so interpreting the operators ˆ a † and ˆ a as creating and annihilating photons of certainmomentum was leapt at very easily – thus their eponymous current names. Soon after, thelanguage introduced by the quantization of the electromagnetic field crept into all physicalproblems where the Hamiltonian could be recast into the likeness of Eq. (31). In everyinstance (with the exception of the harmonic oscillator), the particles which the creation13nd annihilation operators were meant to create or annihilate were given a specific name.In the case of lattice vibrations, for example, the particles are called phonons[9, 10]. In thisnarrative, the operator ˆ a † ˆ a is interpreted as counting the number of particles occupying acertain state, hence its name occupation number operator , usually indicated as ˆ n . In thisview, the state | n = 0 i contains no particle, so we call it the vacuum state . The absence ofparticles in the vacuum state though does not necessarily imply a state of zero energy. Asit mostly happens, it is characterized by a specific energy. For the harmonic oscillator, Eq.(15), this vacuum energy corresponds to E = ~ ω . And anytime a particle is added to thesystem, the energy of the latter increases by ∆ E = ~ ω – which is commonly referred to as a quanta of energy . In addition, we can imagine creating any state | n i from the vacuum state | i . It is easy to prove from Eq. (30) that: | n i = (cid:0) ˆ a † (cid:1) n √ n ! | i . (32)The vacuum state of the harmonic oscillator can thus accommodate any finite number ofparticles, hence the particles must be bosons. The normalization factor √ n ! can be inter-preted as accounting for the indistinguishability of these bosons. (It is interesting to notethe analogy here with the solution to the Gibbs paradox.). The particle states {| n i} are theorthonormal basis elements of a vector space called the Fock space [9–11]. In this space, twobasis states | n i and | n ′ i are orthogonal to each other if they differ in their occupation num-bers. The Fock space for the harmonic oscillator is bosonic because, as stated above, we canfill the vacuum state with any number of bosons. There are Fock spaces where one cannotfill the vacuum state with more than one particle. These are called fermionic Fock spaces,and the particles in question are fermions. Finally, we mention that from the properties ofthe operators ˆ a † and ˆ a , Eq. (30), one derives the commutation relation: (cid:2) ˆ a, ˆ a † (cid:3) = ˆ I (33)which can also be easily derived applying the definition of ˆ a † and ˆ a given at Eq. (29).In discussing the harmonic oscillator above, we made a very deep conceptual leap when wetransitioned from first to second quantization. This has to do with the interpretation of theinteger n . Under the first quantization scheme, n just indexed eigenfunctions like ψ n ( x ) andtheir respective energies (Eq. (16)). But according to the second quantization scheme, wecame to see the same integer n as being an eigenvalue of the operator ˆ a † ˆ a and also indicates14he number of bosons occupying a Fock space ket. Both interpretations are correct and canbe used interchangeably. However, some caution is needed in how far we drag the meaning of n together with the creation and annihilation operators in second quantization. In general,the particles created or annihilated according to second quantization represent excitationsin the system under study. This is quite clear from the relation between n and the quantumharmonic oscillator’s eigenenergy, Eq. (15): n indicates how energetically excited the stateof the oscillator is. And the notion of ˆ a † and ˆ a creating and annihilating some bosons,respectively, according to the second quantization scheme, is just a mathematical constructwhich provide an alternative way of talking about these same harmonic excitations. Thesame applies to lattice vibrations: phonons, like the bosons for the harmonic oscillator, arejust mathematically constructed particles which occupation numbers represent excitationsin lattice vibrations.However, there are also many instances whereby these particles one get from secondquantizing a system are not just the fruits of some mathetmatical trickery, but are realparticles to reckon with in Nature. For example, photons represent the excitations in theelectromagnetic field according to second quantization, and are real. The Higgs bosonrepresent excitations in the Higgs field, and it has recently been detected experimentally.The electromagnetic field and the Higgs field are examples of what we call quantum fields .In fact, all the elementary particles in physics (including the electron) are excitations of aparticular quantum field. The study of quantum fields and their excitations is the subjectof quantum field theory (QFT)[10]. We are not going to need QFT in the discussions thatfollow, but it is important we bear in mind the episteme related to these particles whichtranspire through second quantization.To conclude this brief discussion on the harmonic oscillator and second quantization, weconsider a collection of noninteracting harmonic oscillators, say N in total. If we work inthe occupation number representation, the Hamiltonian here is a direct generalization of Eq.(25): ˆ H = N X i =1 (cid:18) ˆ a † i ˆ a i + 12 (cid:19) ~ ω i (34)It is then very easy to see that a generic state of the system can be expressed as | n , n , . . . , n N i where n i indicates the number of bosons present in the i − th harmonic oscillator. The (over-all) vacuum state of the system is the state in which n i = 0 for any given oscillator i , i.e.15 , , . . . , i , which we shall simply indicate as | i . The eigenenergies are also easily foundto be: E n ,n ,...,n N = N X i =1 (cid:18) n i + 12 (cid:19) ~ ω i . (35)The operators ˆ a † i and ˆ a i operate only on the Fock space of the i − th oscillator, therefore suchoperators with different indexes commute: [ˆ a i , ˆ a i ′ ] = h ˆ a † i , ˆ a † i ′ i = ˆ0 , h ˆ a i , ˆ a † i ′ i = δ i,i ′ ˆ I . (36)As in the case of a single harmonic oscillator, the generic state | n , n , . . . , n N i can begenerated from the vacuum state | i : | n , n , . . . , n N i = N Y i =1 (ˆ a † i ) n i √ n i ! | i . (37)The generic multi-harmonic oscillator state | n , n , . . . , n N i is simply the direct product | n i ⊗ | n i ⊗ . . . ⊗ | n N i . The Fock space for this collection of oscillators is thus the vectorspace tensor product of the Fock spaces of the separated oscillators. The collection of integers { n , n , . . . , n N } is an ordered multiset: the first element denotes the number of bosons inthe first oscillator, the second element is the number of bosons in the second oscillator andso on. Each individual oscillator has its own vacuum state, and each of these has its ownenergy content. The energy ( E ) of the overall vacuum state, | i , is the sum of the energyof the various local vacuum states: E = P Ni =1 ~ ω i . B. The Holstein-Primakoff Transformation
Having discussed the harmonic oscillator, we are now ready to discuss the Holstein-Primakoff (HP) transformation.Given an arbitrary particle of spin- j ( j = , , , , , . . . ), we commonly represent its spinquantum states through a set of orthonormal states {| j, m i} which are each a simultaneouseigenstate of the operators ˆ J ˆ J ˆ J and ˆ J z : ˆ J ˆ J ˆ J | j, m i = j ( j + 1) ~ | j, m i (38a) ˆ J z | j, m i = m ~ | j, m i (38b)As we know, the possible values of the magnetic spin quantum number, m , depends on j and the nature of j . If j is an integral integer, then m = ± j, ± ( j − , . . . , , while for16alf-integral j , m = ± j, ± ( j − , . . . , ± / . The fact that the values of m can be either: 1)all integral or, 2) all half-integral integers, creates some discomfort when writing algorithmsfor doing computations on systems which may potentially involve spin quantum numbers ofdifferent types. The ideal way to go about this would be to have a simple way of representingthe possible m values as function of some parameter which is independent of whether j isintegral or half-integral. In fact, for a given j , it is easy to verify that the possible values of m can be simply expressed as: m = j − n , where n = 0 , , , . . . , j . (39)Note that, alternatively, we could also have chosen the form: m = j + n ′ , where n ′ = 0 , − , − , . . . , − j . (40)However, the form given in (39) is preferred because the parameter n admits only nonnegativeintegers. But another reason why we choose n (Eq. 39) over n ′ (Eq. 40) is that the formerallows a very natural transition to the second quantization scheme, and, as we shall seeshortly, this is what leads to the Holstein-Primakoff transformation.Indeed, the relation (39) implies that there is one-to-one correspondence between m and n for fixed j . So instead of the eigenstates {| j, m i} , we can equally make use of the states {| j, n i} – which are also simultaneous eigenstates of ˆ J ˆ J ˆ J and ˆ J z . Just as the states {| j, m i} are often simply indicated as {| m i} , we will often indicate the states {| j, n i} as {| n i} . Forexample, for a spin- / , the states | m = +1 / i and | m = − / i in the | m i representationbecome | n = 0 i and | n = 1 i in the | n i representation, respectively. If the new states {| n i} bring to mind the simple quantum harmonic oscillator, then you already get where this isgoing. In particular, if we combine Eq. (38) and Eq. (39), we get: ˆ J z | n i = ~ ( j − n ) | n i (41)where we have effected the transformation | j, m i → | n i . Repeating the reasoning which ledto Eq. (25), we see that we may express ˆ J z as ˆ J z = ~ (cid:16) j − ˆ b † ˆ b (cid:17) (42)where ˆ b † ˆ b | n i = n | n i . (43)17n the language of occupation number representation, the operator ˆ b † ˆ b is the occupationnumber operator for some particles. The vacuum state of these particles, i.e. | n = 0 i , isseen to correspond to the state | j, j i ; and as we increase n , m decreases by the same degree.To delineate further the parallelism between this way of talking about spin states andthe harmonic oscillator, let us consider a single spin − j interacting with a static magneticfield B o = B o e z . The spin Hamiltonian, as we know, is: ˆ H = − γB o ˆ J z (44)where γ is the spin’s gyromagnetic ratio. According to Eq. (42), we can rewrite thisHamiltonian also as: ˆ H = (cid:16) j − ˆ b † ˆ b (cid:17) ~ ω (45)where ω := − γB o is the Larmor frequency. And the eigenvalues of ˆ H are easily seen to be: E n = ( j − n ) ~ ω . (46)The quanta of energy is still ~ ω and the vacuum energy here is thus E = j ~ ω . We alsoobserve from Eq. (46) that the occupation number for the particles which ˆ b † ˆ b counts indicatesexcitations in the spin, just as we saw for the quantum harmonic oscillator (Eq. (15)). Itis worthwhile to point out that while the vacuum state of the harmonic oscillator is a trueground state (that is, the state with the lowest energy), in the present case, the vacuumstate | i is not always the ground state. Here, | i is the ground state only when γ > ,while for γ < the | i corresponds to the highest excited state.Eqs. (46) and (45) rightly reminds us of Eqs. (31) and (15) from the harmonic oscillatorproblem, respectively. In addition, the operators ˆ b † and ˆ b have the same properties as the ˆ a † and ˆ a encountered when we discussed the harmonic oscillator, respectively, Eq. (30);namely: ˆ b | n i = √ n | n − i , ≤ n ≤ j (47a) ˆ b † | n i = √ n + 1 | n + 1 i , ≤ n ≤ j − (47b)where the ranges on n have been set so as to remain consistent with Eq. (39). The abovebounds imposed on n mark a very important difference between the new states {| n i} for thespin and those for the quantum harmonic oscillator. We shall come back to this point verysoon. 18q. (42) gives a second quantization representation of the spin operator ˆ J z , but this isnot enough to allow us to fully do spin dynamics in second quantization. We also need toexpress the operators ˆ J x and ˆ J y in terms of the operators ˆ b † and ˆ b . To achieve this, it israther convenient to deal with their linear combinations ˆ J ± = ˆ J x ± i ˆ J y . From the theoryof angular momentum, we know that, for example, ˆ J + | m i = ~ p ( j − m )( j + m + 1) | m + 1 i . (48)From Eq. (39), we have that if | m i → | n i , then | m + 1 i → | n − i , therefore Eq. (48) inthe occupation number representation becomes: ˆ J + | n i ≡ f + (ˆ b † , ˆ b ) | n i = ~ p j − ( n − √ n | n − i (49)where f + (ˆ b † , ˆ b ) is simply the operator ˆ J + written in terms of ˆ b † and ˆ b . Our objective is tofind f + (ˆ b † , ˆ b ) . From the first equation of Eq. (47a), we note that: f + (ˆ b † , ˆ b ) | n i = ~ p j − ( n −
1) ˆ b | n i . (50)Given that ˆ b | n i ∝ | n − i , and the final state must remain | n − i , the operator which gen-erates the coefficient p j − ( n − must be in function of the occupation number operator, ˆ b † ˆ b , Eq. (43). We are then lead to the conclusion that: ˆ J + = ~ q j − ˆ b † ˆ b ˆ b . (51)Since ˆ J − and ˆ J + are Hermitian conjugate of each other, it follows immediately from Eq.(51) that: ˆ J − = ~ ˆ b † q j − ˆ b † ˆ b . (52)Eqs. (42), (51) and (52) constitute the Holstein-Primakoff transformation . In the HPtransformation the usual spin operators ˆ J x , ˆ J y , ˆ J z are all written in function of a singleoperator ˆ b † and its Hermitian conjugate, ˆ b . The particles which ˆ b † ( ˆ b ) creates (annihilates)are called the Holstein-Primakoff bosons . We emphasize that the HP bosons simply representspin excitations, and we should not go beyond this interpretation. Interestingly, while ˆ J + | m i ∝ | m + 1 i and ˆ J − | m i ∝ | m − i , we observe from Eqs. (51) and (52) that ˆ J + | n i ∝| n − i and ˆ J − | n i ∝ | n + 1 i , which is consistent with the fact that an increase in m impliesa decrease in n and vice versa (because the sum m + n is conserved, Eq. (39)). Therefore, the19perator ˆ J + (unlike ˆ b † ) annihilates HP bosons, but ˆ J − creates them. In complete analogyto the m − representation where we know ˆ J + cannot increase the magnetic spin quantumnumber m of the state | m i indefinitely, in the occupation number representation ˆ J + cannotannihilate the HP bosons indefinitely but ends when n = 0 , which corresponds to the HPvacuum state. Analogously, ˆ J − can fill the vacuum state with a maximum number of HPbosons, namely, n = 2 j . This limit is set by the operator q j − ˆ b † ˆ b in the definition of ˆ J − ,Eq. (52). This is all consistent with Eq. (39) where it was evident that n is a nonnegativeinteger, and whose range is bounded by the spin quantum number j : ≤ n ≤ j .From Eqs. (42), (51) and (52) we also derive the following: ˆ J x = ~ q j − ˆ b † ˆ b ˆ b + ˆ b † q j − ˆ b † ˆ b (53a) ˆ J y = ~ q j − ˆ b † ˆ b ˆ b − ˆ b † q j − ˆ b † ˆ b i (53b) ˆ J ˆ J ˆ J = ~ j ( j + 1)ˆ I . (53c)Most importantly, one can verify that the Holstein-Primakoff transformation preserves thecommutation relations: h ˆ J α , ˆ J β i = i ~ ǫ αβγ ˆ J γ (54)or their equivalent: h ˆ J z , ˆ J ± i = ± ~ ˆ J ± (55a) h ˆ J + , ˆ J − i = 2 ~ ˆ J z (55b)where α, β and γ represent any of the three directions x, y, z and ǫ αβγ is the three-dimensionalLevi-Civita symbol. Hence, working in the HP representation is just the same as in the m − representation, in the sense that the physics does not change.In analogy to the quantum harmonic oscillator, the HP spin states {| n i} can all begenerated from the vacuum state | i : | n i = (ˆ b † ) n √ n ! | i , ≤ n ≤ j . (56)We could also obtain | n i from | i by applying repeatedly ˆ J − . This has the advantage ofincorporating inevitably the bound on n . Certainly, | n i ∝ ˆ J − ~ ! n | i . (57)20ndeed, the following identities can be easily proved: | n i = (cid:16) ˆ J − ~ (cid:17) n n ! q(cid:0) jn (cid:1) | i (58) | n i = Q nk =1 (2 j + k − ˆ b † ˆ b ) / n ! q(cid:0) jn (cid:1) (ˆ b † ) n | i (59) | n i = ˆΛ( j, n ) (ˆ b † ) n √ n ! | i (60)where, the operator ˆΛ( j, n ) is defined as: ˆΛ( j, n ) := r(cid:16)(cid:16) j +1 − ˆ b † ˆ bn (cid:17)(cid:17)q(cid:0) jn (cid:1) (61)and where (cid:0)(cid:0) xn (cid:1)(cid:1) is defined as: (cid:16)(cid:16) xn (cid:17)(cid:17) := x ( x + 1)( x + 2) . . . ( x + n − n ! . (62) x in Eq. (62) can be a number or an operator. For a nonnegatve integer x , (cid:0)(cid:0) xn (cid:1)(cid:1) is whatwe call the multiset coefficient (read as " x multichoose n "), and it gives the number ofcombinations of length n we can get from a set of x elements, if we allow repetition ofelements and disregard order[6]. For example, the multisets of cardinality we can get froma set of cardinality like { i, j, k } are: { i, i } , { i, j } , { i, k } , { j, j } , { j, k } , { k, k } ; so (cid:0)(cid:0) (cid:1)(cid:1) = 6 ,which can be checked from Eq. (62). By definition, (cid:0)(cid:0) xn (cid:1)(cid:1) = 1 when n = 0 and x is a scalar;while (cid:0)(cid:0) xn (cid:1)(cid:1) = ˆ I when n = 0 and x is an operator. To prove Eqs. (59) and (60) from (58) thefollowing identity is very useful: (ˆ b † ) k q j − ˆ b † ˆ b = q j + k − ˆ b † ˆ b (ˆ b † ) k (63)for nonnegative integer k . The identity can be easily proved using the properties of theoperators ˆ b † and ˆ b , Eqs. (47a).The operator ˆΛ( j, n ) defined in Eq. (61) is a very interesting one. It can be easily shownthat: ˆΛ( j, n ) | n i = | n i , if ≤ n ≤ j nonexistent , otherwise . (64)21q. (60) is very telling: (ˆ b † ) n √ n ! | i certainly yields the quantum harmonic oscillator state | n i – Eq. (32). So, for now, n can be any nonnegative integer, just like with the quantumharmonic oscillator. But by multiplying (ˆ b † ) n √ n ! | i from the left by the operator ˆΛ( j, n ) , werestrict the range of values for n to ≤ n ≤ j . Evidently, ˆΛ( j, n ) plays the role of a"cut-off" operator: while it forbids n from being outside the range ≤ n ≤ j , it reducesto the identity operator when n is within the same range. Interestingly, this observation onthe properties of ˆΛ( j, n ) entails a very peculiar implication. It begins with the observationthat, lim j →∞ ˆΛ( j, n ) = ˆ I (65)according to Eq. (61). But if ˆΛ( j, n ) becomes the identity operator for any occupationnumber n , then there is a one-to-one correspondence between the states {| n i} in Eq. (60)and the quantum harmonic oscillator states {| n i} we saw in Eq. (32). This means that theeigenvectors of the quantum harmonic oscillator are the same as those of a spin- ∞ particlein Fock space[3]. Conversely, we can interpret the Fock space in which the spin states {| n i} are defined as a truncated version of the quantum harmonic’s Fock space by means of thethe operator ˆΛ( j, n ) .The next sections dedicated to multispin systems notwithstanding, we briefly introducethe topic here. Given an arbitrary finite collection of spins, their spin quantum numbersform a multiset, which we shall indicate as A , i.e. A = { j , j , . . . , j N } . Multispin statesare normally indicated in the m − representation either using the coupled or uncoupledrepresentations[12]. In the latter, the states, {| j m , j m , . . . , j N m N i} , are simultaneouseigenstates of the individual ˆ J ˆ J ˆ J i and ˆ J zi operators. In the HP occupation number represen-tation, | j m , j m , . . . , j N m N i 7→ | n , n , . . . , n N i , where n i , ≤ n i ≤ j i , is the numberof HP bosons introduced into the HP vacuum state of the i − th spin, and N is the totalnumber of spins. Each possible combination of occupation numbers { n , n , . . . , n N } alsoforms an ordered multiset and corresponds to a specific multispin state. The overall HPvacuum state is obviously the | n , n , . . . , n N i state with n i = 0 for every i . If we indicatesuch state as | i , it trivially follows from Eq. (60) that: | n , n , . . . , n N i = ˆΛ( j , n ; . . . ; j N , n N ) N Y i =1 (ˆ b † i ) n i √ n i ! | i (66)22here, ˆΛ( j , n ; . . . ; j N , n N ) := N Y i =1 r(cid:16)(cid:16) j i +1 − ˆ b † i ˆ b i n i (cid:17)(cid:17)q(cid:0) j i n i (cid:1) (67)If we consider a multiset of noninteracting spins A = { j , j , . . . , j N } placed in a regionwith static magnetic field B o = B o e z , the multispin Hamiltonian of the system is a facilegeneralization of Eq. (45): ˆ H = N X i =1 (cid:16) j i − ˆ b † i ˆ b i (cid:17) ~ ω i (68)and the eigenenergies are simply: E n ,n ,...,n N = N X i =1 ( j i − n i ) ~ ω i . (69)and the total vacuum energy is E = P i j i ~ ω i . Like the single spin HP vacuum state, themultispin analogue is not necessarily the ground state. To get the ground state from | i ,one needs to fill all localized vacuum states characterized by nonnegative ω i (i.e. negative γ i ) with HP bosons to their maximum capacity and leave those with negative ω i empty. IV. MULTISPIN MAGNETIC RESONANCE WITH THE HOLSTEIN-PRIMAKOFFTRANSFORMATION
We have seen that the uncoupled representation of the multispin states associated with A = { j , j , . . . , j N } are {| n , n , . . . , n N i} in the HP occupation number representation, andthe possible occupation numbers ( n , n , . . . , n N ) are multisets of integers. To effectivelydiscuss multispin dynamics, it is convenient to have a reasonable shorthand for indicating ageneric state | n , n , . . . , n N i . The so-called index compression maps are devised just for thispurpose and they will be our first order of business in this section. In the second part, wewill analyze the eigendecomposition problem for some general forms of isotropic multispinHamiltonians using the HP representation. A. Index compression maps
A detailed discussion on index compression maps is beyond the scope of this article, sowe are going to keep the mathematical details to the minimum. We will concern ourselves23ere with only one special kind of index compression maps, denoted η . For more on indexcompression maps see [13]. Given a spin multiset A = { j , j , . . . , j N } , the index compressionmap, η , maps each possible combination of occupation numbers ( n , n , . . . , n N ) to a uniqueinteger of the set { n } . The mapping is one-to-one. η is simply a function which has theHP occupation numbers n , n , . . . , n N as variables. To be more specific: n = η ( n , n , . . . , n N ) = N X i =1 W R,i · n i . (70)where, W R,i := δ N,i + (1 − δ N,i ) N − i Y k =1 d i + k (71)where δ i,i ′ is the Kronecker delta and d i is the dimension of the spin- i ’s spin Hilbert space,i.e. d i = 2 j i + 1 . In more comprehensible terms, for i = N , W R,i is the product of all the d i ′ with i ′ > i ; while for i = N , W R,i = 1 . For instance, W R, = d · d · · · d N . In thesame spirit, W R, then corresponds to the dimension D H of the entire spin Hilbert spaceof the multispin system. We thus infer from Eq. (70) that the possible values of n are: , , , . . . , ( D H − [13]. Note that since each of the states | n , n , . . . , n N i constitute a basisfor the multispin Hilbert space, the dimension of the latter must coincide with the number ofall possible combinations of the HP occupation numbers { n , n , . . . , n N } , which is precisely D H . These combinations are ordered and restricted: restricted in the sense that for each n i , ≤ n i ≤ j i , as we mentioned above. The range of n is thus consistent with the dimensionof the multispin Hilbert space.We adopt the following notation: if n is related to the multiset { n , n , . . . , n N } throughthe index compression map η , Eq. (70), then we shall write the ket | n , n , . . . , n N i as | n i .The integers n will be indicated with the font O , , , , . . . to differentiate them from thesingle spin states. The complete HP vacuum state of the system always has n = , and so | , , . . . , i 7→ | i , consistent with our previous notation for this particular state.Let us return to the deuterated hydroxymethyl radical ( • CH OD). There are only fournonzero spins we need to take into consideration for the simulation of its ESR spectrum: 1)the radical electron, 2) the first hydrogen nuclear, 3) the second hydrogen nuclear, and 4)the deuterium nuclear. If we number the spins according to the same order, then the spinmultiset in this case is A = { j , j , j , j } = { , , , } . And from Eq (70) we can write: n = 12 n + 6 n + 3 n + n . (72)24n table I the elements of • CH OD’s spin Hilbert space basis in the uncoupled represen-tation ( | m , m , m , m i ) are translated into the HP occupation number representation( | n , n , n , n i ), and their respective shorthand notation | n i according to the index compres-sion map η is also given. As one can easily observe, the relationship between the integers n and the occupation numbers ( n , n , n , n ) according to the mapping η strictly dependson how we label the spins.It is interesting to note that once we know the labels on the spins from the multiset A = { j , j , . . . , j N } , we can determine the corresponding HP occupation numbers n , n , . . . , n N for any given n . This is because the map η is invertible and has a unique inverse: η − ( n ) = { n , n , . . . , n N } . In this sense, Eq. (70) can be viewed as a linear multivariable Diophantineequation in the variables n , n , . . . , n N . The problem is somehow complicated by the factthat each variable n i is restricted to a specific range. In any event, there is a simple algorithmone can follow to solve Eq. (70) for n , n , . . . , n N when n is fixed. It goes as follows: Findthe integer n , ≤ n ≤ j , which maximizes the product W R, n but still keeps the latterless or equal to the given n ; find then n in the range ≤ n ≤ j which maximizes the sum W R, · n + W R, · n , still keeping the sum less or equal to n . Continue with the procedure tillyou get to n N . In the particular case where along the way you get to a certain occupationnumber n k , where n ≤ n k < n N , such that n = W R, · n + W R, · n + . . . + W R,k · n k , then n k +1 = n k +2 = . . . = n N = 0 .For instance, say we want to find the HP occupation numbers n , n , n , n correspondingto n = for the • CH OD radical (assuming we maintain the same choice of spin labellingas done above). This amounts to solving the Diophantine equation: = 12 n + 6 n + 3 n + n (73)in accordance with Eq. (72), knowing that n , n , n can be or or , while n can beany nonnegative integer not greater than . According to the algorithm discussed above,we see that n must necessarily be , n and n must be and, finally, n = 2 . Thus, | i = | , , , i , which is in agreement with Table (I). It therefore goes without sayingthat by means of the map η , we can easily encode multispin states in the form of integers,and the decoding process is as easy as that of encoding.Before moving on, we point out that, unlike the multispin HP boson Fock states, the map η cannot be applied to the boson Fock states for a finite collection of quantum harmonic25 m , m , m , m i | n , n , n , n i | n i (cid:12)(cid:12) + , + , + , +1 (cid:11) | , , , i | i (cid:12)(cid:12) + , + , + , (cid:11) | , , , i | i (cid:12)(cid:12) + , + , + , − (cid:11) | , , , i | i (cid:12)(cid:12) + , + , − , +1 (cid:11) | , , , i | i (cid:12)(cid:12) + , + , − , (cid:11) | , , , i | i (cid:12)(cid:12) + , + , − , − (cid:11) | , , , i | i (cid:12)(cid:12) + , − , + , +1 (cid:11) | , , , i | i (cid:12)(cid:12) + , − , + , (cid:11) | , , , i | i (cid:12)(cid:12) + , − , + , − (cid:11) | , , , i | i (cid:12)(cid:12) + , − , − , +1 (cid:11) | , , , i | i (cid:12)(cid:12) + , − , − , (cid:11) | , , , i | i (cid:12)(cid:12) + , − , − , − (cid:11) | , , , i | i (cid:12)(cid:12) − , + , + , +1 (cid:11) | , , , i | i (cid:12)(cid:12) − , + , + , (cid:11) | , , , i | i (cid:12)(cid:12) − , + , + , − (cid:11) | , , , i | i (cid:12)(cid:12) − , + , − , +1 (cid:11) | , , , i | i (cid:12)(cid:12) − , + , − , (cid:11) | , , , i | i (cid:12)(cid:12) − , + , − , − (cid:11) | , , , i | i (cid:12)(cid:12) − , − , + , +1 (cid:11) | , , , i | i (cid:12)(cid:12) − , − , + , (cid:11) | , , , i | i (cid:12)(cid:12) − , − , + , − (cid:11) | , , , i | i (cid:12)(cid:12) − , − , − , +1 (cid:11) | , , , i | i (cid:12)(cid:12) − , − , − , (cid:11) | , , , i | i (cid:12)(cid:12) − , − , − , − (cid:11) | , , , i | i Table I: Basis kets of • CH OD’s spin Hilbert space in the uncoupled representation( | m , m , m , m i ), the HP representation ( | n , n , n , n i ) and their shorthand notation | n i according to the index compression map η .oscillators (not even for a single oscillator). This is because the single spin HP boson Fock26paces we encounter in Nature are always finite in dimension (due to the action of theoperator ˆΛ( j, n ) , Eq. (64)) – this is in stark contrast to the boson Fock space related to asingle quantum harmonic oscillator, which is of infinite dimension. B. Eigendecomposition of Isotropic Multispin Hamiltonians
We discuss here how the HP transformation and the index compression map η can bedeployed to solve the eigendecomposition problem for isotropic multispin Hamiltonians. Wewill be led to a computational scheme that fails to scale polynomially with the number ofspins in some cases, but yet far better than the tout court approach. One can considerwhat we discuss below as the first step to systematically address the curse of dimensionalitywithin the occupation number representation approach, with particular emphasis on the HPtransformation. It is our opinion that numerical analysis of the problem alone won’t get usfar enough and that we should resort to the latter after we have exhausted our mathematicaltricks and conveniences within a given theoretical approach. One of the advantages of the HPtransformation is that it has the tendency of making it possible to find analytical solutionsto seemingly intractable problems. This analytical element can be very useful when writingalgorithms for multispin systems.Symmetry arguments cannot be ignored in problems like the one we are about to discuss,and we will greatly make use of them. But it is important to note that symmetry argumentscan get us to confidently reach certain conclusions without sharing light on how to effectivelycarry on the needed calculations. This is where the HP transformation and the indexcompression map η come to our aid. To make the symmetry arguments more intelligiblewe are not going to consider the general isotropic multispin Hamiltonian in the presence ofa static magnetic field straightaway, but we will get there through a gradual process. Andalong the way, we will discuss various symmetries and analyze the spin Hamiltonians weencounter with the help of the HP transformation and η .Without loss of generality, we have set ~ = 1 in the following.27 . The isotropic multispin Hamiltonian in the absence of an external magnetic field In the introduction section of the current paper, we saw that in the absence of ex-ternal magnetic fields, the isotropic spin Hamiltonian of an arbitrary collection of spins A = { j , j , . . . , j N } takes the form of ˆ H spin − spin , Eq. (1), which we now rewrite as[1]: ˆ H spin − spin = X i>i ′ T i,i ′ ˆ J ˆ J ˆ J i · ˆ J ˆ J ˆ J i ′ = X i>i ′ T i,i ′ (cid:20) ˆ J zi ˆ J zi ′ + 12 (cid:16) ˆ J + i ˆ J − i ′ + ˆ J − i ˆ J + i ′ (cid:17)(cid:21) . (74) H spin − spin manifestly remains the same in all directions (rotationally invariant). Not sur-prisingly, it is a zero rank spherical tensor operator. Another way of proving this is byrecalling the definition of spherical tensor operators. A spherical tensor operator of rank k can be defined as an operator with (2 k + 1) irreducible components which obey the followingcommutation relations[9]: h ˆ J z , ˆ T ( k ) q i = ~ q ˆ T ( k ) q (75a) h ˆ J ± , ˆ T ( k ) q i = ~ p ( k ∓ q )( k ± q + 1) ˆ T ( k ) q ± (75b)with k = 0 , , , . . . and − k ≤ q ≤ + k . (Spherical tensor operators are often definedaccording to how they behave under rotation[12]. This involves the use of Wigner matricesand we would like to avoid them for now. We therefore stick to the definition above, whichis much simpler to handle.) For example, if we compare the commutation relations in Eqs.(55) with Eqs. (75), we note that the trio ˆ J z , ˆ J ± must be proportional to the componentsof a rank k = 1 spherical tensor. Indeed, if we set: ˆ T (1)0 = ˆ J z ˆ T (1) ± = ∓ √ J ± (76)we see that the relations in Eq. (75) are perfectly observed.Going back to ˆ H spin − spin , it can easily be verified that: h ˆ J ztot , ˆ H spin − spin i = ˆ0 (77a) h ˆ J ± tot , ˆ H spin − spin i = ˆ0 (77b)where, ˆ J ztot = P Ni =1 ˆ J zi and ˆ J ± tot = P Ni =1 ˆ J ± i . In line with Eqs. (75), we observe that since ˆ H spin − spin is not a null operator, Eqs. (77) hold only if ˆ H spin − spin is proportional to a zerorank tensor ( k = 0 ). From Eq. (77b), we also derive that: h ˆ J xtot , ˆ H spin − spin i = ˆ0 (78a) h ˆ J ytot , ˆ H spin − spin i = ˆ0 (78b)28he conclusion we draw from the commutation relations stated in Eqs. (77a), (78a) and(78b) is that the total spin angular momentum must be conserved along the three axes, justas one would expect from a rotationally invariant operator like ˆ H spin − spin . The operators ˆ J xtot , ˆ J ytot , ˆ J ztot are thus constants of motion for an isolated spin system whose spin Hamil-tonian is given by ˆ H spin − spin . Note that any linear combination of these three operators isalso a constant of motion. Thus, the total spin angular momentum vector operator ˆ J ˆ J ˆ J tot , forexample, is also a constant of motion. Not only that: any operator of the form ˆ X r , where r = 1 , , , . . . and ˆ X is any of the operators ˆ J ˆ J ˆ J tot , ˆ J xtot , ˆ J ytot , ˆ J ztot , is also a constant of mo-tion. The same applies to any linear combination of such powers of operators. However, notall pairs of operators of this vast family of constants of motions commute with each other.For example, ˆ J xtot does not commute with ˆ J ytot . This is crucial because the easiest way todetermine the energy spectrum of ˆ H spin − spin , for example, is to express it as a function of asubset of mutually commuting operators which belong to this infinitely numerable family ofconstants of motion. Unfortunately, as far as we know at the moment, ˆ H spin − spin as given inEq. (74) fails to succumb to this mathematical contrivance due to the arbitrary differencebetween the coupling constants { T i,i ′ } . Numerical diagonalization of H spin − spin , therefore,seems the most reasonable route to the eigenvalues and eigenvectors of ˆ H spin − spin , especiallywhen dealing with multispin systems.As mentioned earlier, one could employ the tout court approach whereby one simplydiagonalizes H spin − spin in whole. The other approach is to first divide H spin − spin into in-dependent sub-units and then proceed with the diagonalization of each sub-unit. The HPtransformation and the index compression map η o discussed in §III B and §IV A, respectively,enables us to break H spin − spin into sub-units in a computationally efficient manner.First of all, consider two arbitrary multispin states | n i (= | n , . . . , n N i ) and | n ′ i (= | n ′ , . . . , n ′ N i ) . From Eq. (77a), we have that: D n (cid:12)(cid:12)(cid:12) h ˆ J ztot , ˆ H spin − spin i (cid:12)(cid:12)(cid:12) n ′ E = 0 . (79)Given that, ˆ J ztot = N X i =1 (cid:16) j i − ˆ b † i ˆ b i (cid:17) = J − ˆ N (80)(where J ≡ P Ni =1 j i is the total spin of the system, and ˆ N ≡ P Ni =1 ˆ b † i ˆ b i the total HPbosons occupation number operator), we note that – as one would expect – the kets {| n i} ˆ J ztot : ˆ J ztot | n , . . . , n N i = ˆ J ztot | n i = ( J − n ) | n i (81)where n = P Ni =1 n i is the total number of HP bosons contained in | n i . In light of Eq. (81),Eq. (79) reduces to the form: ( n − n ′ ) D n (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n ′ E = 0 (82)where it becomes evident that D n (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n ′ E = 0 when n = n ′ . In other words, anecessary but not sufficient condition for the matrix element D n (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n ′ E to be nonzerois that the two multispin states | n i and | n ′ i are represented by the same total number of HPbosons. The total number of HP bosons is thus conserved for a spin system with Hamiltonian ˆ H spin − spin , Eq. (74); i.e. h ˆ N , ˆ H spin − spin i = 0 (this commutation relation easily follows fromEqs. (77a) and (81)).The importance of the relation in Eq. (82) even goes further than what we have con-cluded so far. Indeed, if D n (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n ′ E is identically zero for any pair of states | n i and | n ′ i with different number of total HP bosons, it implies that in the matrix representation of ˆ H spin − spin according to the basis {| n i} , the matrix elements between states with the sametotal number n of HP bosons constitute a block which is orthogonal to all other possibleblocks characterized by a different total number of HP bosons. In other words, if we decom-pose the Hamiltonian ˆ H spin − spin in the basis {| n i} , all the multispin states | n i occupied bythe same total number n of HP bosons compose a subspace B n of the HP bosons Fock space,and subspaces characterized by different n will be orthogonal to each other, i.e. B n ⊥ B n ′ if n = n ′ . But since there is a one-to-one correspondence between the states {| n i} and thenormal multispin states {| j m , . . . , j N m N i} , it means that the normal multispin Hilbertspace of the system is also decomposed into orthogonal subspaces in similar fashion. Inmore concise mathematical terms, the subspace B n is simply the set of kets defined as: B n := {| n i | ˆ N | n i = n | n i , ≤ n ≤ ( D H − } (83)(recall D H is the dimension of the Hilbert space).It is only natural at this point to ask ourselves to determine the range of n for a givenmultiset A of spins. Since n = P Ni =1 n i , and ≤ n i ≤ j i , it is clear that ≤ n ≤ J .30onsequently, for a given collection of spins whose Hamiltonian operator is the q = 0 -th component of a rank k spherical tensor (see Eq. (75)) like ˆ H spin − spin , the system’sHilbert space can always be decomposed into (2 J +1) orthogonal subspaces. This particularobservation is far from new: as a matter of fact, it is well known that whenever an operatorcommutes with ˆ J ztot , the total spin magnetic quantum number M z is conserved – which leadsto the creation of orthogonal subspaces each characterized by a particular M z . Since thetotal spin is J , the range of M z is − J ≤ M z ≤ J , which means there are (2 J + 1) distinctpossible values of M z , hence (2 J + 1) orthogonal subspaces. Indeed, as it follows from Eq.(81), the connection between M z and n is very simple: M z = J − n . Nonetheless, the useof n is far more convenient computationally than M z because, unlike M z whose values canbe either all integers or all half-integers, the possible values of n are always nonnegativeintegers, independent of the collection of spins at hand. This is not the only advantage ofusing the HP transformation, and we shall discuss others shortly.The realization of the orthogonal subspaces {B n } also implies that the matrix represen-tation of ˆ H spin − spin , H spin − spin , is block-diagonalized in the basis {| n i} . In fact, the relation H spin − spin = J M n =0 B n = diag ( B , B , . . . , B J ) = B . . . B J (84)holds true. The block matrix B n collects the matrix elements of ˆ H spin − spin between ketsbelonging to the subspace B n , i.e. B n := nD n (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n ′ E (cid:12)(cid:12)(cid:12) ∀ n , n ′ ∈ B n o . (85)Given a multiset A of spins, we shall indicate the dimension of the block matrix B n as Ω A ,n ,which is also the dimension of the orthogonal subspace B n . Thus, for consistency, D H = J X n =0 Ω A ,n . (86)The very crucial implication of Eq. (84) is that the eigendecomposition of H spin − spin can be equally achieved by eigendecomposing each B n independently, with the benefit of31educing the computational cost (in both computing time and memory). If the compu-tational cost of the eigendecomposition of a square matrix of dimension M is O ( M p ) (inmost cases p ≈ ), we see that the computational cost of eigendecomposing H spin − spin toutcourt is O ( D p H ) . In contrast, if we choose to determine the eigenvalues and eigenvectors of H spin − spin by eigendecomposing the block matrices B n , then the computational cost turnsout to be: O (Ω p A , ) + O (Ω p A , ) + . . . + O (Ω p A , J ) . Due to the relation in Eq. (86), it is patentlyobvious that for p > , O ( D p H ) > O (Ω p A , ) + O (Ω p A , ) + . . . + O (Ω p A , J ) (87)which confirms our assertion that eigendecomposing the block matrices B n costs less thanthe tout court eigendecomposition of H spin − spin . (In this analysis, we have tacitly assumedthat the computational cost of generating the H spin − spin matrix is the same as the total costof generating the block matrices { B n } .) a. Integer partitions. Dimension of the submatrices B n . Density and sparse-ness of H spin − spin Another significant feature of the HP transformation is that it allows us to easily deter-mine the dimension Ω A ,n of the various submatrices B n , analytically. This is a feat hardlyachievable using the normal spin representation.The analytical determination of the value of Ω A ,n has been extensively covered in [3] sowe shall just limit ourselves to some key points in the following. The interested Reader maysee [3] for a more elaborate exposition of the problem.Recall that | n i = | n , n , . . . , n N i , and n = η o ( n , n , . . . , n N ) (see Eq. (70)). As alreadyseen above, if ˆ N | n i = n | n i , then n = P Ni =1 n i (see Eqs. (80) and (81)). We also need tobear in mind that each n i is restricted to the range ≤ n i ≤ j i . Since Ω A ,n is the numberof kets | n i which contain the same total number of HP bosons (Eqs. (83) and (85)), theproblem at hand is equivalent to asking: in how many distinct ways can we distribute n indistinguishable objects among N sites, knowing that the i − th site can contain at most j i objects? The number of ways this can be done is exactly Ω A ,n . Counting problems likethis is the subject of a branch of discrete mathematics called Enumerative Combinatorics[6],which fundamentally deals with how to count the elements of a finite set. Several ways ofcomputing Ω A ,n are illustrated in [3] but we shall focus here on just one of these, namely32he generating function approach. Generating functions in Enumerative Combinatorics are formal power series (in one or multiple variables) whose coefficients are proportional to thesolutions to a counting problem.Consider for example the partition of integers. The partition of a positive integer n isa way one can obtain n through the sum of positive integers, order being irrelevant. Forexample, take the integer , since: (88)the sums on the RHS are all partitions of . Each summand in a given partition is referredto as a part . Thus, the partitions and have both two parts, while the partition has three parts. Say p ( n ) the total number of partitions of the integer n . p ( n ) is known as the partition function . We see that p ( n = 4) = 5 , for example, since theinteger can be partitioned in five different ways as illustrated above. We may then posethe following problem: given an arbitrary (positive) integer n , in how many ways can wepartition it? Or, in other words, what is the value of its p ( n ) ? Variant forms of thisproblem have been proposed throughout centuries. According to known records, it seemsGottfried Wilhelm von Leibniz (1646-1716) was the first to pose a variant of this problem ina letter to Johann Bernoulli (1667-1748)[14]. Leibniz was interested in how many ways aninteger could be partitioned into two, three, etc., parts. In Leibniz’s problem, we clearly seethere is a constrain on the number of parts. This falls under what we call today restrictedpartitions . For example, there is only one way to partition into three parts, and that is . In a letter to the great Leonhard Euler, one Naudé wanted to know how manyways the integer 50 can be obtained from the sum of seven parts which are unequal to eachother[15–17]. Here, the restriction is both on the nature of the parts (i.e. each part cannotappear more than once) and the total number of parts. It is not far fetched to assume thatNaudé’s letter is what led the great Euler to his memorable Observationes analyticae variaede combinationibus [ Various analytical observations about combinations ] presented to theSt. Petersburg Academy in 1741 but published in 1751[15]. This is where Euler introduced33or the first time the concept of generating function (though he did not coin the name) inthe theory of partitions (and by extension, number theory). For example, he showed that p ( n ) is the coefficient of q n when Q ∞ i =1 11 − q i is expanded in powers of q . Namely[16, 17], ∞ Y i =1 − q i = ∞ X n =0 p ( n ) q n = 1 + q + 2 q + 3 q + 5 q + 7 q + 11 q + 15 q + 22 q + . . . (89)The function Q ∞ i =1 11 − q i is thus said to be the generating function for p ( n ) . Recall at thebeginning of this section we defined generating functions as being formal power series. Thisis so because their variables have no intrinsic meaning. In the case of Eq. (89), q is thevariable but though its exponents and coefficients have definite interpretations in relation toour counting problem, q is devoid of any. Perhaps, the most famous generating function inNumber theory and Combinatorics is the one reported in Eq. (89). In response to Naudé’sproblem, Euler showed that if we denote with ˜ p k ( n ) the number of ways n can be partitionedinto k mutually unequal parts, then[15, 16]: q k ( k +1) / k Y i =1 − q i = ∞ X n =0 ˜ p k ( n ) q n . (90)Thus, q k ( k +1) / Q ki =1 11 − q i is the generating function for ˜ p k ( n ) . To solve Naudé’s originalproblem, we set n = 50 and k = 7 ; so from Eq. (90) it follows that ˜ p (50) = 522 , which isobtained by expanding q Q i =1 11 − q i and taking the coefficient of the term q . Hence, thereare ways one can write the integer as the sum of exactly seven mutually unequalpositive integers. Euler’s original paper, now translated into English by Jordan Bell[16],is an excellent and easy-to-read introduction to generating functions and it is accessible toanyone with a high school knowledge of algebra. We strongly recommend it to Readers whoare new to the concept.Instead of determining the solution to counting problems by means of generating func-tions, we may seek explicit formulae for the desired quantities. For instance, we may askourselves if there is an explicit formula for p ( n ) . It is a problem which quite a number ofgenerations of mathematicians had to wrestle with. The first significant breakthrough camein Godfrey H. Hardy (1877-1947) and Srinivasa Ramanujan’s (1887-1920) celebrated 1918paper entitled Asymptotic formulae in combinatory analysis [18]. In that paper, Hardy andRamanujan presented an asymptotic series for p ( n ) but whose main defect was that it failed34o converge. A couple of decades had to pass before Hans Rademacher (1892-1969) derivedan improved version of Hardy and Ramanujan’s result which, finally, had no convergenceproblem[19].The Reader might be wondering why we have dedicated so many lines to integer parti-tions. The reason is very simple: the determination of the value of Ω A ,n may be reinterpretedas an integer partition problem. This is how the close relation between multispin dynam-ics, on one hand, and Number theory and Enumerative Combinatorics (specifically in thiscase, the theory of partitions) elegantly emerges from the HP transformation. Indeed, if, asnoted earlier, Ω A ,n is the number of ways the integer n can be obtained through the sum n = P Ni =1 n i , where ≤ n i ≤ j i , it is clear that: Ω A ,n is the number of partitions of the integer n into N parts, with the i − thpart restricted to the range ≤ n i ≤ j i .Certainly, just like ˜ p k ( n ) , Ω A ,n is restricted to a fixed number of parts. But unlike p ( n ) and ˜ p k ( n ) , Ω A ,n has the following properties:1. the minimum value of a part in any of its partitions is , not ;2. the partitions of interest here are all ordered (so for example, must beconsidered different from ); this is because each part n i refers to a distinctspin.Consider for example the deuterated hydroxymethyl radical ( • CH OD), which correspondsto the multiset of spins A = { j , j , j , j } = { , , , } (see §IV A). Therefore, according tothe conventions discussed in §IV A, we know that since j = j = j = and ≤ n i ≤ j i ,then ≤ n , n , n ≤ , while ≤ n ≤ . Ω A ,n is then the number of ordered partitions of n into exactly N = 4 parts, such that the first three parts are at most and the fourth partis at most (with being an admissible value of a part). We also need to bear in mind that Ω A ,n is also equivalent to the number of states | n i which contain exactly n HP bosons. Forexample, with n = 2 , we find from table I that the states | n i which contain in total two HPbosons are: | i = | , , , i | i = | , , , i | i = | , , , i| i = | , , , i | i = | , , , i | i = | , , , i| i = | , , , i • CH OD, Ω A , = 7 . That is, there are seven ways of partitioningthe integer into four parts, with the restriction that the first three parts cannot exceedthe value of 1 and the fourth part cannot be greater than 2. Thus, the subspace B n =2 for • CH OD is of dimension Ω A , = 7 . In the particular case of • CH OD, n = 2 corresponds tothe total spin magnetic number M z = J − n = − . Note that these results actuallystill hold for any multispin system with three spin- and one spin- .How can we determine the value of Ω A ,n without having to explicitly write down all thekets | n i (and how the HP bosons are disposed in them)? This problem has been solved anddiscussed in [3]. We discuss here only the generating function for Ω A ,n without going intomuch details. The interested Reader may see [3] for further discussions, explicit formulaefor Ω A ,n and proofs.Before we present the generating function for Ω A ,n it is advisable we briefly see whatis meant by the q-analogue of a nonnegative integer n . The q − analogue of the integer n ,indicated as [ n ] q , is defined as[3]: [ n ] q := 1 − q n − q = 1 + q + q + . . . + q n − . (91)Note that lim q → [ n ] q = n .It can be shown that the generating function for Ω A ,n , G A , Ω ( q ) , is[3]: G A , Ω ( q ) = Y α (cid:16) [2 j α + 1] q (cid:17) N α (92)where the index α runs over distinct values of the spin multiset A and N α is the multiplicityof the α − th distinct element in A . Hence, Y α (cid:16) [2 j α + 1] q (cid:17) N α = J X n =0 Ω A ,n q n . (93)It is worth noting that if we take the limit q → of Eq. (93) we get Eq. (86). Toillustrate the use of Eq. (93), let us consider once again • CH OD. We know that for thisradical, A = { , , , } . A has thus only two distinct elements: j α =1 = and j α =2 = 1 .The multiplicity of j α =1 is while that of j α =2 is . Therefore, N α =1 = 3 and N α =2 = 1 .Moreover, [2 j α =1 + 1] q = [2] q = (1 + q ) and [2 j α =2 + 1] q = [3] q = (1 + q + q ) . Thus, in thecase of the radical • CH OD, Eq. (93) becomes: (1 + q ) (1 + q + q ) = X n =0 Ω A ,n q n . (94)36ince, (1 + q ) (1 + q + q ) = 1 + 4 q + 7 q + 7 q + 4 q + q (95)we conclude that: Ω A , = Ω A , = 1 Ω A , = Ω A , = 4 Ω A , = Ω A , = 7 (96)which is in agreement with the previous value we found for Ω A , . In regards to • CH OD, itssubspaces B n , their respective dimension and basis elements are reported in table II. Whilethe correspondence between the first two columns will always remain the same, the basiskets | n i spanning each subspace depends on the order chosen when labelling the spins. Thekets in table II are the same kets in table I, therefore they correspond to the same orderedmultiset of spins A = { j , j , j , j } = { , , , } . Subspace B n Dimension of subspace, Ω A ,n Basis elements, | n iB | iB | i , | i , | i , | iB | i , | i , | i , | i , | i , | i , | iB | i , | i , | i , | i , | i , | i , | iB | i , | i , | i , | iB | i Table II: Orthogonal subspaces B n , their respective dimension and basis kets for theradical • CH OD in the case whereby the system’s spin Hamiltonian is proportional to the q = 0 − th component of a spherical tensor of rank k – ˆ H spin − spin , Eq. (74), being anexample.What immediately catches our attention looking at the RHS of Eq. (95) (or the secondcolumn of table II) is the "palindromic" distribution of Ω A ,n ’s values. This is not an exceptionbut the rule. Indeed, one can prove that the generating function G A , Ω ( q ) is a reciprocal polynomial. Namely[3]: Ω A ,n = Ω A , J − n . (97)A polynomial P ( x ) = a + a x + a x + . . . + a s x s with real coefficients is said to be reciprocal(or palindromic ) if its coefficients are such that a i = a s − i for all i = 0 , , . . . , s [20]. An equally37quivalent definition is that P ( x ) is reciprocal if x s P (cid:0) x (cid:1) = P ( x ) [20]. The implication of(97) is that the subspaces B n and B J − n have the same dimension. We shall give a simpleproof of Eq. (97) later when we discuss time-reversal symmetry. For now, we also point outthat G A , Ω ( q ) is also unimodal [20], namely, there exists an integer n ′ such that: Ω A , ≤ Ω A , ≤ Ω A , ≤ . . . ≤ Ω A ,n ′ ≥ Ω A ,n ′ +1 ≥ Ω A ,n ′ +2 ≥ . . . ≥ Ω A , J . (98)This n ′ is found to be n ′ = ⌊ J ⌋ [3], where ⌊•⌋ is the floor function. Clearly, Ω A , ⌊ J ⌋ is themaximum value the { Ω A ,n } may assume. Note that Ω A , ⌊ J ⌋ may not be the only Ω A ,n withthis maximum value (afterall, G A , Ω ( q ) is reciprocal). For example, in the case of • CH OD,whose A = {
12 3 , } , J = 5 / ; hence, ⌊ J ⌋ = ⌊ ⌋ = 2 . Therefore, we expect Ω A , to have themaximum value any Ω A ,n here can possibly have. In fact, that is the case but Ω A , is not theonly Ω A ,n with this maximum value since Ω A , also has the same value as Ω A , (Eq. (96)).For more on the properties of G A , Ω ( q ) see [3].In the limit case whereby all the spins are spin- , i.e. the spin system is univariate (see§II) with j = , A = n N o , the generating function G A , Ω ( q ) takes the simple form: G A , Ω ( q ) = (1 + q ) N = N X n =0 Ω A ,n q n (99)Consequently, the dimension of the subspace B n , Ω A ,n = Ω n N o ,n , in this limit case is givenby a binomial coefficient: Ω n N o ,n = (cid:18) Nn (cid:19) . (100)Let us consider for instance a spin system which consists of N = 10 spin- s, i.e. A = n
12 10 o .We are dealing here with a univariate spin system. If the spin Hamiltonian of the systemcommutes with ˆ J ztot (so the Hamiltonian is proportional to a spherical tensor of rank k but q = 0 , see Eq. (75)), the tout court approach will have us diagonalize a matrix of dimension = 1024 . But with the creation of the subspaces B n , we know that G A , Ω ( q ) = (1 + q ) = 1 + 10 q + 45 q + 120 q + 210 q + 252 q + 210 q + 120 q + 45 q + 10 q + q (101)38hus, instead of eigendecomposing a matrix of dimension to determine the eigenvec-tors and eigenvalues of the spin Hamiltonian, we can equally obtain the same results byeigendecomposing smaller matrices of dimension , , , , and .The distribution of the values of Ω A ,n does not only inform us about the dimension ofthe subspaces B n , but they also help us to quantify how sparse the matrix H spin − spin is. Letus define the density ζ ( A ) of an arbitrary matrix A as the ratio between the number of itsnonzero elements and the dimension of A . Analogously, we define the sparseness χ ( A ) of thematrix A as the fraction of the elements of the latter which are identically zero. Naturally,for any given matrix A : ζ ( A ) + χ ( A ) = 1 . (102)It is easy to prove that ζ ( H spin − spin ) is subject to the tight upper bound: ζ ( H spin − spin ) ≤ P J n =0 Ω A ,n D H (103)independent of the specific values of the coupling constants T i,i ′ . If we go back to the radical • CH OD, for example, ζ ( H spin − spin ) ≤ ∼ . . This means no matter what the values ofthe constants T i,i ′ are, the matrix representation of the Hamiltonian ˆ H spin − spin cannot havemore than of its elements being nonzero. This is not only true for • CH OD, but alsofor all multiset of spins A = n
12 3 , o . Indeed, Eq. (103) holds for any arbitrary multiset ofspins whose spin Hamiltonian is proportional to the zero-th component of a spherical tensorof rank k , where k = 0 , , , . . . For a univariate spin system of N spin- , i.e. A = n N o , it follows from Eqs. (100) and(103) that: ζ ( H spin − spin ) ≤ P Nn =0 (cid:0) Nn (cid:1) N = (cid:0) NN (cid:1) N (104)from which one derives that for very large N , ζ ( H spin − spin ) . √ πN . (105)Eq. (105) shows in unambiguous terms that for a system like A = n N o , ζ ( H spin − spin ) tendsto zero as N becomes very large. For instance, if we consider the spin system A = n
12 1000 o ,i.e. a collection of spins, all of spin- , then ζ ( H spin − spin ) . √ π ≈ . ; whichmeans the elements of H spin − spin which are identically zero will never drop below ofthe total entries. 39nterestingly, we note that if we begin with a spin system A = n N o , and substituteone of the spins with a particle whose spin quantum number is greater than / , the new H spin − spin is less denser than the original. We can thus imagine creating any collection of N spins from the multiset A = n N o by substitution. Given that anytime we substitutea spin- / with a greater spin the density ζ ( H spin − spin ) reduces, it implies that for anyimaginable multiset A = n N o of spins, ζ ( H spin − spin ) < (cid:0) NN (cid:1) N (106)and in the limit of a large number of spins, ζ ( H spin − spin ) < N − δ , δ = 12 (1 + log N π ) . (107) b. Eigenenergies of a system of equivalent spins in the absence of an externalfield. We illustrate here another powerful application of the HP transformation. In the limit casewhereby the multiset A = { j , j , . . . , j N } consists of solely N interacting equivalent spins(§II), ˆ H spin − spin (Eq. (74)) reduces to the form: ˆ H spin − spin = T X i>i ′ ˆ J ˆ J ˆ J i · ˆ J ˆ J ˆ J i ′ = T X i>i ′ (cid:20) ˆ J zi ˆ J zi ′ + 12 (cid:16) ˆ J + i ˆ J − i ′ + ˆ J − i ˆ J + i ′ (cid:17)(cid:21) . (108)Since, ˆ J ˆ J ˆ J tot = X i ˆ J ˆ J ˆ J i + 2 X i>i ′ ˆ J ˆ J ˆ J i · ˆ J ˆ J ˆ J i ′ , (109) ˆ H spin − spin in Eq. (108) may be written as: ˆ H spin − spin = T ˆ J ˆ J ˆ J tot − X i ˆ J ˆ J ˆ J i ! . (110)Given that ˆ J ˆ J ˆ J tot commutes with ˆ J ˆ J ˆ J i , we can easily determine the eigenvalues of ˆ H spin − spin ifwe know those of ˆ J ˆ J ˆ J tot and ˆ J ˆ J ˆ J i . Fortunately, the expression for the eigenvalues of this set ofcommuting operators is well known. From Eq. (53c), for example, we have that: ˆ H spin − spin = T " J tot ( J tot + 1) − X i j i ( j i + 1) ˆ I . (111)As simple as Eq. (111) may seem, the actual computation of the eigenenergies is not aneasy task for a generic A with N > . The difficulty here lies in computing the various40otal spin angular momentum J tot and their multiplicities. But the J tot and their respectivemultiplicities are the Clebsch-Gordan series, which is the multiset of all the possible total(spin) angular momenta one can get by coupling the N elements of A . For example, if N = 2 (in which case Eqs. (108) and (111) hold irrespective of whether the two spins areequivalent or not), we well know that | j − j | ≤ J tot ≤ j + j . Certainly, when N > one could determine the Clebsch-Gordan series by coupling j and j , and then couple theresulting angular momenta with j , followed by the coupling of the new set of resultingangular momenta with j , and one repeats the scheme till one gets to j N . Needless to say,this is truly cumbersome. It is understandable that if one has an easy and computationallyefficient method of computing the Clebsch-Gordan series for a generic collection of spins, onealso reduces dramatically the computational cost of determining the eigenvalues of ˆ H spin − spin in Eq. (108). Here again, the HP transformation is invaluable.The problem of determining analytically the Clebsch-Gordan series for an arbitrary col-lection of spins A has been solved in [3]. We will therefore limit ourselves here to citing somesalient results from [3] without providing any proof.Consider an arbitrary collection of spins A = { j , j , . . . , j N } . The addition of the spinangular momentum of the elements of A will result in a collection (or multiset) of total spinangular momenta { J tot } . Let us indicate the distinct elements of { J tot } as J , J , J , . . . , J m ,with J being the maximum, J the second highest and so forth. Consecutive distinct J tot differ by , hence: J κ = J − κ , where κ = 0 , , , . . . , m . (112)Clearly, J = j + j + . . . + j N . The first big challenge we encounter here is computing thevalue of J m . It can be proved that[3]: J m = υ A · H ( υ A ) + (1 − H ( υ A )) · (2 J mod 2)2 (113)where, υ A := 2 · max A − J (114)and where H ( x ) is the Heaviside step function, defined here to be H ( x ) := , if x < , if x ≥ . (115)41ay N J the total number of distinct elements in { J tot } . Evidently, N J = J − J m + 1 . Atthis point, we may represent the multiset { J tot } as: { J tot } = (cid:8) J λ , J λ , . . . , J λ m m (cid:9) , where λ κ is the multiplicity of J κ .Having determined all the distinct elements of { J tot } thanks to Eqs. (112) and (113), weare half-way through in getting to the Clebsch-Gordan series. All we need to do now is todetermine the multiplicities { λ κ } . Once again, we can easily solve the problem by makinguse of the HP transformation. One can show that the generating function for λ κ , G A ,λ ( q ) ,is related to G A , Ω ( q ) – Eq. (92) – through the relation[3]: G A ,λ ( q ) = (1 − q ) G A , Ω ( q ) = J +1 X κ =0 λ κ q κ . (116)The polynomial G A ,λ ( q ) is antipalindromic since, λ κ = − λ J +1 − κ . (117)Moreover, it can be easily proved that: m X κ =0 λ κ = Ω A ,m (118)i.e. the cardinality of the Clebsch-Gordan series is exactly Ω A ,m .In our current quest to determine the eigenvalues of the Hamiltonian ˆ H spin − spin of Eq.(108), we see that since the multiset A is fixed, the eigenvalues of ˆ H spin − spin can be distin-guished on the basis of J tot . We may thus indicate these eigenenergies as E J tot , or E κ on thebasis of Eq. (112). We shall employ the latter in the following. It then follows from Eqs.(111) and (112) that: E κ = E − T κ (2 J + 1 − κ ) (119)where κ = 0 , , , . . . , m = ( J − J m ) , and E := T J − X i j i ! . (120)The following observations on the system readily follows from Eq. (119):1. if T > , then E m is the ground state ( J tot = J m ) energy and E is the energy of thehighest excited state; 42. if T < , then E is the energy of the ground state ( J tot = J ) , and E m is the energyof the highest excited state;3. E κ +1 − E κ = − T J κ . This means that the energy of multispin states with consecutive J tot at the lower end of { J tot } (i.e. as J κ approaches J m ) are relatively less spacedcompared to their counterparts at the higher end.Furthermore, given that the degeneracy of J κ is (2 J κ + 1) , the total degeneracy of the energylevel E κ will thus be given by the product: λ κ (2 J κ + 1) . The sum total of the degeneracy ofthe energy levels must return the dimension of the system’s spin Hilbert space, D H . Thisleads us to the relation: D H = m X κ =0 λ κ (2 J κ + 1) (121)which is a simple sum rule one can deploy to spot flaws in the calculations. Other sum rulescan be derived from Eq. (121) by combining it with Eqs. (112) and (118).To illustrate the usefulness and potential of the relations and techniques discussed above,let us consider the multiset of spins A = n
12 7 , o . Assuming A is a multiset of equivalentspins whose Hamiltonian is given by (111), we ask: what are the energy levels of the systemand their respective degeneracy? As simple as this problem may appear, it is rather difficult –if not computationally time consuming – to solve using the conventional spin representation.The dimension of the spin Hilbert space alone is D H = 3456 , though less than . of thematrix elements of the spin Hamiltonian is nonzero. It would be a waste of resources toconstruct such a matrix and then eigendecompose it to find the energy levels. It is evenmore computationally challenging to directly employ Eq. (111) to compute the energy levelsbecause – to the best of our knowledge – there is not, hitherto, a generally valid efficientalgorithm to compute the Clebsch-Gordan series for arbitrary multispin systems. In theframework of the HP transformation, such a problem can be effortlessly solved without theneed to eigendecompose H spin − spin .Naturally, the eigenenergies { E κ } of the system is given by Eq. (119). Concerningthe distinct values of the multiset { J tot } , we easily determine its maximum to be J = . From Eq. (113), the minimum J tot is found to be J m = . Thus, there are N J = J − J m + 1 = 7 distinct values of J tot , and m = J − J m = 6 . The distinct J tot are: { J , J , J , J , J , J , J } = (cid:8) , , , , , , (cid:9) . We now determine the multiplicity of each J tot by employing the generating function G A ,λ ( q ) , Eq. (116). First of all, from Eq. (92), we43ave that G A , Ω ( q ) = (1 + q ) (1 + q + q ) . Therefore, G A ,λ ( q ) = (1 − q )(1 + q ) (1 + q + q ) = 1 + 9 q + 38 q + 99 q + 174 q + 207 q + 145 q + . . . − q = X κ =0 λ κ q κ (122)We therefore conclude that the Clebsch-Gordan series for A = n
12 7 , o is given by themultiset { J tot } = n ,
112 9 ,
92 38 ,
72 99 ,
52 174 ,
32 207 ,
12 145 o . With these values, we can now easilycompute all the eigenenergies E κ of the system from Eq. (119). The results are reported intable III, where we have also calculated the degeneracy of each energy level. As expected onthe basis of Eq. (121), the sum of these degeneracies gives exactly – which we recall isthe dimension of the spin Hilbert space of the system. κ J κ λ κ E κ − E Degeneracy of E κ − T − T − T − T − T − T P = 3456 Table III: Energy levels and their respective degeneracy for a system of equivalent spinscomposed of seven spin- / and three spin- , A = n
12 7 , o , whose Hamiltonian is given byEq. (108).It is worth considering the same problem but this time with A = n N o . Surely, J = N/ ,and from Eq. (113), we derive that, J m = (2 J mod 2)2 if N > if N = 2 . N > , we observe that J m = 0 if N is even, while J m = 1 / if N is odd. Thus, thenumber of distinct J tot to expect is N J = ⌈ N − ⌉ + 1 , and m = ⌈ N − ⌉ , where ⌈•⌉ is the ceilingfunction. Therefore, the distinct J tot we get from the Clebsch-Gordan series for A = n N o are: { J , J , . . . , J ⌈ N − ⌉ } = { N/ , N/ − , . . . , N/ − ⌈ N − ⌉} . Here, the multiplicity λ κ of J κ obeys the simple relation: λ κ = (cid:18) Nκ (cid:19) − (cid:18) Nκ − (cid:19) = N + 1 − κN + 1 − κ (cid:18) Nκ (cid:19) , κ = 0 , , , . . . , m . (123)Eq. (123) readily follows from Eqs. (99), (100) and (116). c. Time-reversal symmetry. Creation of submatrices B n . To the best of our knowledge, in the general context of Eq. (74), it is not in general possibleto easily determine the eigenvalues of ˆ H spin − spin like we just did for an arbitrary collection ofequivalent spins in the last section. In general, the eigendecomposition of ˆ H spin − spin requires,first of all, creating the submatrices { B n } and then eigendecomposing each separately, Eq.(84). To create the submatrix B n in the HP representation, we need a general formula for ˆ H spin − spin ’s matrix elements in terms of the HP bosons occupation numbers. Indeed, givenany two basis kets | n ′ i = | n ′ , n ′ , . . . , n ′ N i and | n i = | n , n , . . . , n N i , we derive from Eq. (74)that: D n ′ (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n E = D n ′ , n ′ , . . . , n ′ N (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n , n , . . . , n N E = δ n ′ , n X i>k T i,k ( j i − n i )( j k − n k )+ 12 X i>k T i,k s(cid:18) n i + 1 ± (cid:19) (cid:18) j i − n i + 1 ∓ (cid:19) (cid:18) j k − n k + 1 ± (cid:19) (cid:18) n k + 1 ∓ (cid:19) × ∆ i,k δ n ′ i ,n i ± δ n ′ k ,n k ∓ (124)where, ∆ i,k := Y l = i,k h n ′ l | n l i . (125)In deriving Eq. (124), we have made use of Eqs. (51), (52) and (47a). Eq. (124) confirmsonce again that D n ′ (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n E is identically zero when | n ′ i and | n i differ in the totalnumber of HP bosons they contain.It is worth noting that the creation of the submatrices B n can be greatly simplified anddone in an efficient manner if we make use of the time-reversal symmetry (or T − symmetry).45n the field of magnetic resonance (and in quantum chemistry, in general), we are accus-tomed to space symmetries like rotational and translational invariances, but we seldomlyspeak of time-reversal symmetry. One of the reasons for this, we suppose, can be attributedto the fact that unlike other symmetry operations like rotations and translations which arerepresented by unitary operators in quantum mechanics, T − symmetry is represented byan antiunitary operator (see below) – whose properties depart considerably from the knownunitary operators. In certain areas of quantum chemistry like vibrational spectroscopy, T − symmetry can be of no enlightening use, but in magnetic resonance it is often of vitalimportance as we shall shortly see. For more on T − symmetry, interested Readers may see(in this order): [9, 21, 22]. For an in-depth introduction embedded in some interestingphilosophical discussions, see [23].Let ˆΘ be the time-reversal operator. That is, ˆΘ t ˆΘ − = − t (126)where t is the real parameter which indicates time. Like any operator representing a symme-try transformation, ˆΘ ˆΘ − = ˆΘ − ˆΘ = ˆ . As already remarked above, unlike space rotationor translation operators which are unitary, the time-reversal operator ˆΘ is antiunitary. Themajor difference between these two types of symmetry operators is how they operate oncomplex scalars: unitary operators leave complex scalars intact, while antiunitary oper-ators change complex scalars into their corresponding complex conjugate. For example, ˆΘ e iα ˆΘ − = e − iα ∗ , where α is a complex number whose complex conjugate is α ∗ .There are many interesting properties of the time-reversal operator, but for our purposes,it suffices to know that all angular momentum operators (orbital and spin) are odd underthe operation of time-reversal. Namely[9, 22–24], ˆΘ ˆ J ˆ J ˆ J ˆΘ − = − ˆ J ˆ J ˆ J . (127)This simple relation has many profound consequences which will unfold before us shortly.To begin, it follows from Eq. (127) that ˆΘ ˆ J z ˆΘ − = − ˆ J z . Combining this with Eq. (42),we find that: ˆΘ ˆ b † ˆ b ˆΘ − = 2 j − ˆ b † ˆ b . (128)Till now, we have always seen spin kets of a spin- j as being represented by a number ofHP bosons occupying a certain vacuum space which can accommodate at most j bosons.46nother way of seeing it is to imagine this vacuum state as comprised of j holes (calledHP holes), where each hole can be filled with no more than one HP boson at a time. Inthis picture, the sum of the number of HP bosons and holes is always j . Therefore, if ˆ b † ˆ b counts the number of HP bosons, then (2 j − ˆ b † ˆ b ) counts the number of holes. In other words, (2 j − ˆ b † ˆ b ) is the number operator for the holes. Going back to Eq. (128), the effect of thetime-reversal operator now becomes perspicuous: it transforms the HP occupation numberoperator into a hole number operator. Put in another way, Θ instantly converts HP bosonsinto holes, and vice versa.Before we proceed, let us see how ˆΘ transforms the single spin ket | n i and the multispinket | n i . First of all, we note that Eq. (128) may be rewritten as: h ˆΘ , ˆ b † ˆ b i + = 2 j ˆΘ (129)where h ˆ A, ˆ B i + (:= ˆ A ˆ B + ˆ B ˆ A ) denotes the anticommutation between ˆ A and ˆ B . Say | n i a state ket of a particle of spin- j according to the HP representation, where, as usual, n represents the number of HP bosons. We define the state Θ | n i ≡ | n i as the time-reversed or hole complement of | n i . The state | n i must necessarily be an admissible spin state,else the symmetry operation enacted by Θ would not be reversible. The reversibility of thetime-reversal symmetry operation is guaranteed by the fact that ˆΘ ˆΘ − = ˆΘ − ˆΘ = ˆ I . Aftermultiplying Eq. (129) from the left by | n i , followed by an appropriate rearrangement of theterms, we get: ˆ b † ˆ b | n i = (2 j − n ) | n i (130)where we have made use of Eq. (43). Comparing Eq. (130) with (43), it becomes immedi-ately clear that: ˆΘ | n i ≡ | n i = | j − n i (131)which is, once again, in line with the interpretation that ˆΘ converts HP bosons into holes,and holes into HP bosons. The dual ket relative to ˆΘ | n i is h n | ˆΘ − ≡ h n | .Let us now turn our attention to how ˆΘ operates on a multispin ket | n i . Consider thespin multiset A = { j , j , . . . , j N } . From Eqs. (80) and (128), it follows that: h ˆΘ , ˆ N i + = 2 J ˆΘ (132)which is the analogous of Eq. (129) for multispin systems. Say ˆΘ | n i ≡ | n i the hole com-plement of | n i (= | n , n , . . . , n N i ) . We remind the Reader that the HP bosons occupation47umbers n i in | n , n , . . . , n N i are related to the integer | n i through the index compressionmap η , Eq. (70). If we now multiply Eq. (132) from the right by | n i , we end up with therelation: ˆ N | n i = (2 J − n ) | n i (133)which implies that if | n i belongs to the subspace B n , then its hole complement | n i is anelement of the subspace B J − n . Indeed, since there is a one-to-one correspondence between | n i and | n i , we also conclude that the hole complements of the basis kets in B n constitutethe basis kets of the subspace B J − n . In other words, the two subspaces B n and B J − n are related by time-reversal symmetry, and this holds independent of the nature of the spinHamiltonian. The one-to-one correspondence we just mentioned, implies that the dimensionof the subspaces B n and B J − n must coincide, and indeed, that is something we alreadyknow from Eq. (97).Note that, ˆΘ | n i = ˆΘ | n , n , . . . , n N i = | n , n , . . . , n N i = | j − n , j − n , . . . , j N − n N i . (134)Thus, | n i = | j − n , j − n , . . . , j N − n N i . (135)Hence, by virtue of the index compression map η , it follows from Eq. (70) that: n = N X i =1 W R,i · (2 j i − n i ) = N X i =1 W R,i · (2 j i ) − n . (136)The term P Ni =1 W R,i · (2 j i ) corresponds to the multispin state whereby all the single HPvacuum states are filled with the maximum number of HP bosons they can accommodate.But we know the index compression map η always assigns to this state the integer ( D H − .Thus, n = D H − − n (137)which leads us to conclude that: ˆΘ | n i ≡ | n i = | D H − − n i (138)where we recall once again that D H is the dimension of the multispin Hilbert space. If wetake the • CH OD radical, for example, D H = 24 , so | n i = | − n i . This means that, for48xample, the multispin kets | i and | i (= | ¯ i ) are related by time-reversal symmetry. Infact, from table I, one observes that if one converts all the HP bosons in | i into holes, andvice-versa, one obtains | i . It is worth noting that for a given collection of spins A , Eq.(138) is always valid, independent of how we choose to number-label the spins in A . Thecorresponding dual ket of ˆΘ | n i is h n | ˆΘ − ≡ h n | = h D H − − n | .We are now ready to see what the time-reversal symmetry reveals about the multispinsystem. To understand the importance of time-reversal symmetry here, it is worth consid-ering how it transforms the multispin Hamiltonian ˆ H spin − spin . From Eq. (74) and (127), itfollows that: ˆΘ ˆ H spin − spin ˆΘ − = ˆ H spin − spin (139)(note that the coupling constants T i,i ′ are all real) which means that ˆ H spin − spin is invariantunder time-reversal. This is exactly what we should expect since ˆ H spin − spin describes theHamiltonian of an isolated system, and we know that for such systems the homogeneity oftime applies (i.e. the energy – or, in general, the physics – of the system remains the sameunder any time translation), so ˆ H spin − spin must certainly commute with ˆΘ . Consider nowthe matrix element between the bra h n ′ | and ket | n i according to Eq. (139): D n ′ (cid:12)(cid:12)(cid:12) ˆΘ − ˆ H spin − spin ˆΘ (cid:12)(cid:12)(cid:12) n E = D n ′ (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n E (140a) D n ′ (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n E = D n ′ (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n E . (140b)Eq. (140b) is of vital importance: it tells us that if we know the matrix element D n ′ (cid:12)(cid:12)(cid:12) ˆ H spin − spin (cid:12)(cid:12)(cid:12) n E between the kets | n i and | n ′ i which are elements of the subspace B n , then we also know thecorresponding matrix element between their hole complements in the subspace B J − n . Thisobservation reduces significantly the computational cost of creating the submatrices { B n } :once we know B n we can easily create B J − n . Another interpretation of Eq. (140b) is thatif we convert all HP bosons into HP holes, and vice versa, the physics remain the same.This is a property we may call particle-hole transformation invariance .Eq. (139) has an even more profound implication. Since ˆ H spin − spin is decomposable intosubspaces according to Eq. (84), each eigenvector of ˆ H spin − spin will also belong to onlyone of these subspaces. Say (cid:12)(cid:12)(cid:12) E ( n ) µ E an eigenvector of ˆ H spin − spin but which belongs to thesubspace B n : ˆ H spin − spin (cid:12)(cid:12) E ( n ) µ (cid:11) = E ( n ) µ (cid:12)(cid:12) E ( n ) µ (cid:11) (141)49here E ( n ) µ is the eigenvalue of ˆ H spin − spin relative to (cid:12)(cid:12)(cid:12) E ( n ) µ E . The index µ numbers theeigenvectors of ˆ H spin − spin in B n , so µ = 1 , , . . . , Ω A ,n – where Ω A ,n is the dimension of B n ,§IV B 1 a. As usual, the time-reversed state of (cid:12)(cid:12)(cid:12) E ( n ) µ E , (cid:12)(cid:12)(cid:12) E ( n ) µ E , is obtained by operating ˆΘ on (cid:12)(cid:12)(cid:12) E ( n ) µ E , i.e. (cid:12)(cid:12)(cid:12) E ( n ) µ E = ˆΘ (cid:12)(cid:12)(cid:12) E ( n ) µ E . Given that (cid:12)(cid:12)(cid:12) E ( n ) µ E = (cid:12)(cid:12)(cid:12) E (2 J − n ) µ E , (cid:12)(cid:12)(cid:12) E ( n ) µ E necessarily belongsto the subspace B J − n . Since ˆΘ commutes with ˆ H spin − spin according to Eq. (139), it followsthat: h ˆ H spin − spin , ˆΘ i (cid:12)(cid:12) E ( n ) µ (cid:11) = 0 (142a) (cid:0) E ( n ) µ − E ( n ) µ (cid:1) (cid:12)(cid:12) E ( n ) µ (cid:11) = 0 (142b)from which we deduce that E ( n ) µ = E ( n ) µ . The implication of Eq.(142b) is this: two differenteigenvectors of ˆ H spin − spin related by time-reversal symmetry also share the same eigenvalue.A far more reaching conclusion is that once we are able to create and eigendecompose thesubmatrice B n of H spin − spin , Eq.(84), we do not need to create its hole complement B J − n and eigendecompose it de novo, because the eigenvalues of the two submatrices coincide andtheir eigenvectors are related through ˆΘ . Finding the hole complement of an eigenvector canbe easily done thanks to Eq. (138). For example, let the spin system whose Hamiltonianis given by ˆ H spin − spin in Eq. (74) be the • CH OD radical. From table I, we know that ageneric normalized eigenvector (cid:12)(cid:12)(cid:12) E (1) µ E of the subspace B is given by the linear combination: (cid:12)(cid:12) E ( n =1) µ (cid:11) = c µ, | i + c µ, | i + c µ, | i + c µ, | i (143)where the c µ,i are real coefficients, and P i c µ,i = 1 . The corresponding hole complement of (cid:12)(cid:12)(cid:12) E ( n =1) µ E , (cid:12)(cid:12)(cid:12) E ( n =4) µ E , is an eigenvector of B , and (cid:12)(cid:12) E ( n =4) µ (cid:11) = ˆΘ (cid:12)(cid:12) E ( n =1) µ (cid:11) = c µ, ˆΘ | i + c µ, ˆΘ | i + c µ, ˆΘ | i + c µ, ˆΘ | i = c µ, | i + c µ, | i + c µ, | i + c µ, | i (144)where, in getting to the last step, we applied Eq. (138).It is worth mentioning that there are instances whereby a subspace B n coincides withits time-reversed counterpart. Recall that each subspace B n is characterized by the totalnumber n of HP bosons each of its basis kets contains. So, if each basis ket of the holecomplement of B n , i.e. B J − n , contains (2 J − n ) HP bosons, then B n = B J − n when both50ubspaces are characterized by the same number of HP bosons. That is, when n = 2 J − n ,which happens only when n = J . But n can assume the value of J only when J is aninteger. This means that for multispin systems whose total spin quantum number J is aninteger, the subspace B n = J is its own hole complement. While the eigenvalues of all theother subspaces are atleast double degenerate due to Eq. (142b), the same cannot be saidof B n = J . Conversely, we also conclude that for multispin systems whose total spin quantumnumber J is a half-integer, if the system’s Hamiltonian is invariant under time-reversalsymmetry, then the energy eigenvalues of the system are all atleast double degenerate. Thisconclusion is perfectly inline with Kramer’s degeneracy theorem[9, 24].
2. The isotropic multispin Hamiltonian in the presence of an external static magnetic field
The theoretical machinery developed above are still useful and relevant when we sub-ject the previously isolated multispin system to a static magnetic field B o . The new spinHamiltonian of the system, ˆ H o , is given by Eq. (2). If we take the direction of the exter-nal magnetic field B o to be the axis of quantization e z , so that B o = B o e z , then Eq. (2)becomes: ˆ H o = ˆ Z + ˆ H spin − spin (145)where ˆ H spin − spin is still given by Eq. (74), and ˆ Z := − B o X i γ i ˆ J zi . (146)We know from Eq. (76) that ˆ J z is the zeroth component of a rank k = 1 spherical tensor,so ˆ Z is also the zeroth component of a rank k = 1 tensor. We have also already seen that ˆ H spin − spin is a zero rank spherical tensor. Thus, ˆ H o is the sum of the zeroth components ofspherical tensors of different ranks. But we know from Eq. (75a) that ˆ J ztot commutes withthe zeroth component ( q = 0) of any spherical tensor operator. It therefore follows fromthese considerations that: h ˆ J ztot , ˆ H o i = ˆ0 (147)which is in complete analogy to Eq. (77a). This is significant because it means that all theresults we derived above for the isolated multispin system on the basis of Eq. (77a) alsoapply here. For example, ˆ H o , like ˆ H spin − spin , conserves the total number of HP bosons. So,51 n ′ (cid:12)(cid:12)(cid:12) ˆ H o (cid:12)(cid:12)(cid:12) n E = 0 , if the multispin kets | n ′ i and | n i do not contain the same total number ofHP bosons. The conservation of the total number of HP bosons also implies that, just likein the case of ˆ H spin − spin , ˆ H o subdivides the system’s Hilbert space into (2 J + 1) subspaces: B , B , . . . , B J . Hence, the matrix representation of ˆ H o , H o , can be written in the blockdiagonalized form: H o = J M n =0 B n = diag ( B , B , . . . , B J ) = B . . . B J (148)just as we saw for H spin − spin at Eq. (84), where the definition for B n – Eq. (83) – stillapplies; and B n , in analogy to Eq. (85), is now defined as: B n := nD n (cid:12)(cid:12)(cid:12) ˆ H o (cid:12)(cid:12)(cid:12) n ′ E (cid:12)(cid:12)(cid:12) ∀ n , n ′ ∈ B n o . (149)The generating function for the dimensions { Ω A ,n } of the subspaces {B n } is, therefore,still given by G A , Ω ( q ) , Eq. (92).The expression for the generic matrix element of ˆ H o between | n ′ i = | n ′ , n ′ , . . . , n ′ N i and | n i = | n ′ , n ′ , . . . , n ′ N i , follows directly from Eq. (124) and (146): D n ′ (cid:12)(cid:12)(cid:12) ˆ H o (cid:12)(cid:12)(cid:12) n E = D n ′ , n ′ , . . . , n ′ N (cid:12)(cid:12)(cid:12) ˆ H o (cid:12)(cid:12)(cid:12) n , n , . . . , n N E = δ n ′ , n " − B o X i γ i ( j i − n i ) + X i>k T i,k ( j i − n i )( j k − n k ) + 12 X i>k T i,k s(cid:18) n i + 1 ± (cid:19) (cid:18) j i − n i + 1 ∓ (cid:19) (cid:18) j k − n k + 1 ± (cid:19) (cid:18) n k + 1 ∓ (cid:19) × ∆ i,k δ n ′ i ,n i ± δ n ′ k ,n k ∓ (150)where we can once again observe that ˆ H o conserves the total number of HP bosons. Thesubmatrices B n can be easily created by employing Eq. (150).Here too, the subspaces B n and B J − n are related by time-reversal symmetry. However,the striking difference between the isolated ˆ H spin − spin and ˆ H o is how they are transformed52nder the time-reversal operator ˆΘ . Unlike ˆ H spin − spin , ˆ H o is not time-reversal symmetric.In fact, from Eqs. (127) and (139), we derive that: ˆΘ ˆ H o ˆΘ − = ˆ H o − Z . (151)It must be emphasized that the time-reversal operator ˆΘ we have employed so many timesabove, including Eq. (151), is not a universal time-reversal operator but a local one restrictedto the multispin system. Had ˆΘ been the time-reversal operator acting on the whole Uni-verse, then ˆ H o would certainly commute with ˆΘ since the Universe is an isolated system.Similar results also follow if we consider the multispin system together with the externalmagnetic field B o as one isolated system, and ˆΘ their combined time-reversal operator. Byrestricting ˆΘ only to the multispin system, we have inadvertently partitioned the systeminto two parts: a focus system (the spins) and its environment (the external magnetic field),which are both odd (i.e. they change sign) respect to time-reversal operation. This bipartitepartitioning therefore leads to a broken T − symmetry in ˆ H o .Taking the matrix element of Eq. (151) between the two generic multispin states | n ′ i and | n i , we obtain: D n ′ (cid:12)(cid:12)(cid:12) ˆ H o (cid:12)(cid:12)(cid:12) n E = D n ′ (cid:12)(cid:12)(cid:12) ˆ H o (cid:12)(cid:12)(cid:12) n E − δ n ′ , n Z n , n (152)where, Z n , n := D n (cid:12)(cid:12)(cid:12) ˆ Z (cid:12)(cid:12)(cid:12) n E . If | n ′ i and | n i belong to the subspace B n , then | n ′ i and | n i areelements of B J − n . Thus, Eq. (152) provides an easy way to generate the submatrix B J − n when B n is known.An even more important consequence of Eq. (151) is yet to be unveiled. It has to do withthe relationship between the eigenvalues and eigenvectors of B n and those of B J o − n . Weassume B n is not its own hole complement. Say (cid:12)(cid:12)(cid:12) E ( n ) µ E an eigenvector of ˆ H o in the subspace B n : ˆ H o (cid:12)(cid:12) E ( n ) µ (cid:11) = E ( n ) µ (cid:12)(cid:12) E ( n ) µ (cid:11) , µ = 1 , , . . . , Ω A ,n (153)where E ( n ) µ is the eigenvalue of (cid:12)(cid:12)(cid:12) E ( n ) µ E according to ˆ H o . Surely, the hole complement of (cid:12)(cid:12)(cid:12) E ( n ) µ E , that is: ˆΘ (cid:12)(cid:12) E ( n ) µ (cid:11) ≡ (cid:12)(cid:12) E ( n ) µ (cid:11) = (cid:12)(cid:12) E (2 J − n ) µ (cid:11) (154)is also an eigenvector of ˆ H o , but belongs to the subspace B J − n : ˆ H o (cid:12)(cid:12) E ( n ) µ (cid:11) = E ( n ) µ (cid:12)(cid:12) E ( n ) µ (cid:11) , µ = 1 , , . . . , Ω A ,n (155)53here we recall that Ω A ,n = Ω A , J − n according to Eq. (97). Now, from Eq. (151) we easilyderive the commutation relation between ˆ H o and ˆΘ : h ˆ H o , ˆΘ i = 2 ˆ Z ˆΘ . (156)If we now take the matrix element of Eq. (156) between the eigenvector (cid:12)(cid:12)(cid:12) E ( n ) µ E and its holecomplement (cid:12)(cid:12)(cid:12) E ( n ) µ E , we get: D E ( n ) µ (cid:12)(cid:12)(cid:12) h ˆ H o , ˆΘ i (cid:12)(cid:12)(cid:12) E ( n ) µ E = 2 D E ( n ) µ (cid:12)(cid:12)(cid:12) ˆ Z ˆΘ (cid:12)(cid:12)(cid:12) E ( n ) µ E (157a) E ( n ) µ − E ( n ) µ = 2 D E ( n ) µ (cid:12)(cid:12)(cid:12) ˆ Z (cid:12)(cid:12)(cid:12) E ( n ) µ E (157b)from which follows that: E ( n ) µ = E ( n ) µ − D E ( n ) µ (cid:12)(cid:12)(cid:12) ˆ Z (cid:12)(cid:12)(cid:12) E ( n ) µ E . (158)Eq. (158) asserts that we can easily compute the eigenvalue of the eigenvector (cid:12)(cid:12)(cid:12) E ( n ) µ E fromits hole complement (cid:12)(cid:12)(cid:12) E ( n ) µ E if we already know the latter and its eigenvalue. The significanceof Eq. (158) cannot be stressed enough. For example, suppose we have a univariate spinsystem of spin- , i.e. A = n
12 10 o , whose Hamiltonian in the presence of a static magneticfield is given by ˆ H o , Eq. (145). The total spin quantum number for the system is J = 5 , so,in principle, we can determine the eigenvectors and eigenvalues of ˆ H o by eigendecomposingthe submatrices B , B , B , . . . , B , Eq.(148). But by virtue of Eqs. (154) and (158),there is no need to create and eigendecompose all the eleven submatrices: we only need B , B , B , B , B , B to completely determine the eigenvalues and eigenvectors of ˆ H o .This is an enormous simplification.Most multispin systems of chemical interest present one or more groups of equivalentspins. We can further reduce the computational cost of diagonalizing ˆ H o if we exploit thepresence of these groups in the system. The multiset A of spins can be viewed as the sum of multisets: A = φ ( A ) ] g A g (159)(for example, { , , } ⊎ { , , , } = { , , , } ) where φ ( A ) is the number of groups ofequivalent spins present in A ; the index g runs over the groups of equivalent spins, and A g isthe spin multiset for the g − th group of equivalent spins. Each A g is then transformed into its54lebsch-Gordan series: A g e A g (see §IV B 1 b). Since the elements of the Clebsch-Gordanseries e A g are independent of each other (and irreducible), each element of the Cartesianproduct: φ ( A ) × g e A g = e A × e A × . . . × e A φ ( A ) (160)constitutes an independent multiset of spins subject to the same form of spin Hamiltonianand can be diagonalized independently. Take for example the • CH OD radical. If the ex-perimental conditions are such that the two hydrogen nuclei can be considered as equivalentspins, then the radical consists of three groups of equivalent spins: 1) the unpaired electron, A = (cid:8) (cid:9) ; 2) the two hydrogen nuclei, A = n
12 2 o ; and 3) the deuterium nucleus, A = { } . A and A have only one element each so they coincide with their Clebsch-Gordan multisets e A g . On the other hand, we immediately have that e A = { , } . Hence, e A × e A × e A = (cid:26) (cid:27) × { , } × { } = (cid:26)(cid:26) , , (cid:27) , (cid:26) , (cid:27)(cid:27) . (161)What this means is that when the two hydrogen nuclei in • CH OD are considered equivalent,the radical can be viewed as the sum of two smaller multispin systems, independent of eachother: the first system consists of a spin- / (the unpaired electron) and a spin- particle(the deuterium nucleus), while the second system consists of a spin- / (the same electron)and two spin- particles (the deuterium nucleus and the triplet state of the two hydrogennuclei). The Hamiltonian of both subsystems is still of the form given in Eq. (145). TheHilbert space of both subsystems can be block diagonalized as we saw above. For the firstsubsystem, the dimension of the block matrices are coefficients of the polynomial: (1 + q )(1 + q + q ) = 1 + 2 q + 2 q + q . (162)And the dimension of the block matrices for the second subsystem are the coefficients of thefollowing polynomial: (1 + q )(1 + q + q ) = 1 + 3 q + 5 q + 5 q + 3 q + q . (163)For each subsystem, all that has been discussed above in this section still applies. To reca-pitulate, we see that instead of diagonalizing a matrix of dimension and one of dimension to find the eigenvalues and eigenvectors of ˆ H o for the radical • CH OD as we saw previ-ously, we only need to diagonalize matrices of dimension , and when the two hydrogen55uclei are equivalent. When applied to large spin systems with groups of equivalent spins,this approach reduces the computational cost of diagonalizing the multispin Hamiltoniansignificantly.If we take the naphthalene anion for example, we are dealing with a univariate system of spin-1/2 particles, i.e. A = n
12 9 o . The dimension of the Hilbert space is, therefore, .Suppose we ignore the presence of equivalent spins in the system. Then, to engeindecomposethe systems ˆ H o , we will have to effectively diagonalize four matrices whose dimensionsare , , and , according to the HP transformation scheme discussed above. If theexperimental conditions are sufficiently favorable, the system can be thought of as consistingof three groups of equivalent spins: 1) the electron, A = (cid:8) (cid:9) ; 2) a collection four hydrogennuclei, A = n
12 4 o ; and 3) another collection of four hydrogen nuclei, A = n
12 4 o . Then, e A = { , , } , and e A = { , , } . Thus, e A × e A × e A = (cid:26) (cid:27) × { , , } × { , , } = ((cid:26) , , (cid:27) , (cid:26) , , (cid:27) , (cid:26) , , (cid:27) , (cid:26) , , (cid:27) , (cid:26) , , (cid:27) , (cid:26) , , (cid:27) , (cid:26) , , (cid:27) , (cid:26) , , (cid:27) , (cid:26) , , (cid:27) ) . (164)So the eigendecomposition of the naphthalene anion’s Hamiltonian ˆ H o can be done byconsidering a series of smaller but independent multisets: i) one multispin system comprisedof a spin- / and two spin- particles, i.e. { / , , } ; ii) three multispin systems each of thetype { / , , } , etc. For each of these subsystems, we can apply the HP transformation andcreate the block matrices. In the case of the naphthalene anion, the largest block matrixwe shall encounter comes from the subsystem { / , , } . For this particular subsystem, thedimension of the block matrices follows from the generating function: (1 + q )(1 + q + q + q + q ) = 1 + 3 q + 5 q + 7 q + 9 q + 9 q + 7 q + 5 q + 3 q + q . (165)Hence, if we take into account the two groups of equivalent spins present in the naphtha-lene anion, the largest matrix we will ever have to diagonalize is of dimension , which isremarkable.To fully exploit the presence of groups of equivalent spins in the multispin system, we needthe Clebsch-Gordan coefficients. This constitutes the primary computational challenge in the56ethod just illustrated. One can thus incorporate optimized subroutines for the calculationof Clebsch-Gordan coefficients into ones multispin algorithm. The somehow comfortingobservation we can make here is that, in most of the multispin systems of interest, groups ofequivalent spins hardly exceed six in number – which means a subroutine which can handlethe Clebsch-Gordan coefficients for the addition of up to six angular momenta suffices formost routine computations.
3. Eigendecomposition of Liouvillians
In the normal settings of quantum mechanics, as we have done above, one works in theso-called Hilbert space, where observables like energy, magnetization, etc. are operatorsand the quantum states are vectors. In this way of doing quantum mechanics, the system’senergy are obtained as the expectation values of its Hamiltonian. These energy values are, inprinciple, defined up to a given constant. But in spectroscopic experiments, the propertieswe do observe or measure are essentially related – not to the absolute energies – but tothe difference between these eigenenergies. This fact is nicely conveyed in the resonanceconditions of said experiments. An alternative way of doing quantum mechanics wherebysuch energy differences naturally come up as the expectation values of some operator isthus more suitable to spectroscopy. Such an alternative is accomplished when we formulatequantum mechanics in the so-called Liouville space[25]. In addition, the Liouville spaceformulation makes it relatively easier to treat relaxation processes, compared to the Hilbertspace[25]. It is therefore understandable why the Liouville space formalism is popular invarious fields of spectroscopy, including magnetic resonance. The catch, however, is thatgiven a Hilbert space of finite dimension D H , the dimension of its corresponding Liouvillespace is D H – which roughly translates into even higher computational costs when one worksin Liouville space.As it is well-known, operators in Hilbert space, including the density matrix, become vec-tors ( supervectors ) in Liouville space. Likewise, operations which transformed one operatorinto another in Hilbert space become operators ( superoperators ) in Liouville space.Recall that the Hamiltonian is known to be the generator of the dynamics of the densitymatrix in Hilbert space. When the density matrix becomes a supervector in Liouville space,its dynamics are generated by the superoperator called the Liouvillian . Just as we eigen-57ecompose the Hamiltonian ˆ H in the Hilbert space to get the eigenenergies of the system,the eigendecomposition of its corresponding Liouvillian ˆˆ L returns all possible pairwise dif-ferences between the eigenenergies. The relation between a given Hamiltonian ˆ H and itsLiouvillian ˆˆ L is: ˆˆ L = ˆ H ⊗ ˆ I − ˆ I ⊗ ˆ H ∗ (166)where " ⊗ " denotes the operation of vector space tensor direct product; ˆ I is the identityoperator defined on the same Hilbert space as ˆ H , and ˆ H ∗ is the complex conjugate of ˆ H . We once again observe from Eq. (166) that the dimension of the linear space where ˆˆ L operates, i.e. the Liouville space, is of dimension ( D H ) . For example, if we have a multispinsystem A = n
12 10 o , whose Hamiltonian is given by ˆ H o in Eq. (2), the simulation of themultispin system’s magnetic spectra would require eigendecomposing a matrix of dimension = 1024 according to the tout court approach if we work in Hilbert space. If we chooseto work in the Liouville space, the same tout court approach will have us eigendecomposethe matrix representation of the Liouvillian, which is of dimension = 1 , , .This apparent inconvenience in computational costs one has to grapple with when workingin the Liouville space can be easily overcome. Indeed, if U is the matrix ( supermatrix ) whichdiagonalizes L , i.e. U L U − = D (167)where D is the diagonal supermatrix of the eigenvalues of L , and U is the eigenvector matrixof H : U H U − = D (168)where D is the diagonal matrix containing the eigenvalues of H , then it can be proved that: U = U ⊗ (cid:0) U − (cid:1) T (169)and D = D ⊗ I − I ⊗ D (170)where " ⊗ " in Eqs. (169) and (170) indicates the operation of Kronecker (or matrix direct)product, and I is the identity matrix of the same dimension as D . Eq. (169) is of greatsignificance because it enables us to compute the eigenvector supermatrix U directly from U by means of a simple matrix direct product, without having to diagonalize the Liouvillian L
58e novo. Note that when the Hamiltonian is real, like in the case of the multispin Hamiltonian ˆ H o seen above, U is an orthogonal matrix, and so Eq. (169) reduces to the form: U = U ⊗ U . (171)Combining the relations given in Eq. (169) and (170) with the HP transformation andrelated techniques discussed in previous sections could be very useful in reducing the com-putational cost of simulating the magnetic resonance spectra of multispin systems describedby isotropic Hamiltonians like ˆ H o in Liouville space. For example, going back to our previousexample with the multiset A = n
12 10 o , if we want to work in Liouville space, instead of eigen-decomposing a square matrix of dimension , we can obtain the same eigen-supervectorsand -supervalues through Eq. (171) and (170), respectively, by eigendecomposing onlythe submatrices B , B , B , B , B , B which are of dimension , , , , and ,respectively. The computational cost can be further drastically reduced if we take intoconsideration the presence of groups of equivalent spins as explained in §IV B 2. V. ON SCHWINGER BOSONS
We conclude this paper with a brief introduction to Schwinger bosons, which is anotherkind of spin representation commonly used in condensed matter physics. Schwinger bosonsare closely related to the Holstein-Primakoff bosons, but the two representations are suit-able for certain specific applications. The HP bosons are particularly useful in transformingquantum mechanical problems into counting problems. The Schwinger bosons are extremelypowerful tools when we want to describe the behavior of spin states under unitary transfor-mations like rotations. The literature provides a plethora of applications of the Schwingerbosons, mainly in quantum magnetic studies. This paper won’t be complete without thisbrief introduction to Schwinger bosons. However, the introduction we provide below departsgreatly from what one may find elsewhere. Besides the simplicity of our exposition below,we stress on the close link between Schwinger bosons and HP bosons – which is rarely donein such details in the literature. Readers may see [26, 27] for the mainstream exposition,interpretation and applications of Schwinger bosons.We begin our introduction to Schwinger bosons by noting that the presence of the squareroot of operators in the HP transformation, Eqs. (51) and (52), makes its use in the study59f important problems like the rotation of spin states very inconvenient. The Schwingertransformation, from which derives the Schwinger bosons, provides an alternative way torepresent the spin operators ˆ J ± and ˆ J z in a very simple way without trace of any squareroots. Interestingly, as we show below, the Schwinger representation follows directly fromthe HP transformation.When we introduced the time-reversal operator in §IV B 1 c, Eq. (128), we saw that whenit acts on a spin state represented by HP bosons, it transforms the HP bosons present intoholes and the holes into HP bosons, contemporarily. We also saw that the sum of the numberof holes and HP bosons in any given basis state | n i of a spin- j particle is always j . Andthat, while ˆ b † ˆ b counts the number of HP bosons, the operator (2 j − ˆ b † ˆ b ) counts the numberof holes. The transition from HP transformation to the Schwinger transformation relies onone simple trick: treat the holes and HP bosons as two set of independent particles thatcan be created and annihilated separately (with the caveat that the sum of their occupationnumbers remain constant). If we do so, then we need to assign an occupation numberoperator to the holes. Let us indicate this occupation number operator as ˆ a † ˆ a . Thus, ˆ a † ˆ a = 2 j − ˆ b † ˆ b . (172)Clearly, ˆ a † and ˆ a are the creation and annihilation operators for the holes, respectively; andthey obey the same commutation rules as ˆ b † and ˆ b . Eq. (128) may, therefore, be rewrittenas: ˆΘ ˆ b † ˆ b ˆΘ − = ˆ a † ˆ a . (173)Since { ˆ b, ˆ b † } effect only the HP bosons and { ˆ a, ˆ a † } act on only the holes, the two set ofoperators commute with each other. For the sake of clarity, let ˆ n b ≡ ˆ b † ˆ b and ˆ n a ≡ ˆ a † ˆ a ;thus, the nonnegative integers n b and n a will indicate the number of HP bosons and holes,respectively. We are therefore representing spin states with two types of particles, i.e. HPbosons and (HP) holes:HP representation Schwinger representation | n i 7→ | n b , n a i where, obviously, n = n b . Hence, in the Schwinger representation we need two occupationnumbers to indicate a single spin eigenvector of ˆ J z . Moreover, given that the two sets of60perators { ˆ b, ˆ b † } and { ˆ a, ˆ a † } operate on different particles, it follows that: ˆ b | n b , n a i = √ n b | n b − , n a i ˆ b † | n b , n a i = √ n b + 1 | n b + 1 , n a i (174) ˆ a | n b , n a i = √ n a | n b , n a − i ˆ a † | n b , n a i = √ n a + 1 | n b , n a + 1 i (175)where the nonnegative integers n a and n b are subject to the contraint: n a + n b = 2 j , Eq.(172).In regards to the operator ˆ J z , we note that if we combine Eq. (42) with Eq. (172), weobtain the following expression for ˆ J z in function of ˆ n b and ˆ n a : ˆ J z = 12 (cid:16) ˆ a † ˆ a − ˆ b † ˆ b (cid:17) = 12 (ˆ n a − ˆ n b ) . (176)(recall we have set ~ = 1 ). It is clear from Eq. (176) that the usual spin magnetic number m in | j, m i relates to n a and n b through the expression: m = 12 ( n a − n b ) . (177)Moreover, from Eq. (172), we also have the condition: j = 12 ( n a + n b ) (178)Eqs. (177) and (178), together, constitute the conditions the integers n a and n b must satisfyin order to represent the usual spin state | j, m i .What about the Schwinger representation for ˆ J ± ? To begin, recall the operator ˆ J + = q j − ˆ b † ˆ b ˆ b in the HP representation, Eq. (51). From this expression, we see that theoperator ˆ b first reduces the number of HP bosons by while the operator q j − ˆ b † ˆ b leavesthe resulting state unchanged in the number of HP bosons. The net effect of ˆ J + is, therefore,to reduce n b by . But, given that n a + n b = 2 j , Eq. (178), it follows that if ˆ J + has theeffect of reducing n b by , then it must also have the effect of increasing the number ofholes by the same quantity, i.e. n a n a + 1 . The operator which increases n a by is ˆ a † .Therefore, ˆ J + must be proportional to ˆ a † ˆ b : ˆ J + = c ˆ a † ˆ b = c ˆ b ˆ a † (179)where c is the proportionality constant. If we let ˆ J + operate on | n b , n a i , we find that: ˆ J + | n b , n a i = c ˆ a † ˆ b | n b , n a i (180a) = c p ( n a + 1) √ n b | n b − , n a + 1 i (180b) = c p (2 j − n b + 1) √ n b | n b − , n a + 1 i (180c)61here we have made use of Eqs. (174) and (175). If we compare Eq. (180c) with Eq. (51),we immediately reach the conclusion that c = 1 (bear in mind that n = n b , and as statedabove, we have set ~ = 1 throughout this section). Thus, ˆ J + = ˆ a † ˆ b = ˆ b ˆ a † . (181)Hence, it also follows that: ˆ J − = (cid:16) ˆ J + (cid:17) † = ˆ a ˆ b † = ˆ b † ˆ a . (182a)Eqs. (176), (181) and (182a) constitute the Schwinger transformation. Unlike the HPtransformation, we see that the Schwinger transformation is linear in the operators ˆ a † , ˆ a, ˆ b † , ˆ b .In §III B, we saw that the vacuum state | n = 0 i in the HP representation of the ˆ J z eigen-vectors of spin- j corresponds to the state | j, j i in the normal | j, m i − representation. TheSchwinger representation also has its vacuum state, namely | n b = 0 , n a = 0 i = | , i . Thisvacuum state, unlike the one we encountered in the HP representation, does not correspondto any specific spin state, and its energy is undefined. However, one interesting characteristicof the Schwinger vacuum is that it is the same for all spins. This is in net contrast to theHP vacuum.Just as generic spin states in the HP representation can be created from the HP vacuumstate by filling the latter with a number of HP bosons, generic spin states in the Schwingerrepresentation can be easily created from the Schwinger vacuum | , i by creating a numberof holes and HP bosons. Thus, | j, m i 7→ | n b , n a i = (cid:0) ˆ a † (cid:1) n a √ n a ! (cid:16) ˆ b † (cid:17) n b √ n b ! | , i (183)where the factor √ n a ! n b ! is necessary to keep the state on the RHS normalized like | n b , n a i .If we express n a and n b in Eq. (183) in terms of j and m , using Eqs. (177) and (178), itturns out that: | j, m i 7→ | j − m, j + m i = (cid:0) ˆ a † (cid:1) j + m p ( j + m )! (cid:16) ˆ b † (cid:17) j − m p ( j − m )! | , i . (184)In the Schwinger representation, the HP bosons and HP holes are collectively called Schwinger bosons . 62mong the many useful applications of the Schwinger bosons is the relatively straightfor-ward ease with which they allow the derivation of the famous Majorana formula[28], whichstates the general expression for the transition probability between two generic spin statesof an arbitrary spin in the presence of: 1) a static magnetic field along a specified directionin space (chosen as the quantization axis), and 2) a rf field perpendicular to the static field.Other applications of the Schwinger bosons include the derivation of analytical expressionsfor the matrix elements of the rotation operator. In particular, they can be used to easilyderive a close expression for the elements of the Wigner d − matrix[9, 26]. Closed expressionsfor the Clebsch-Gordan coefficients for the addition of two, three and four angular momentacan also obtained using the Schwinger bosons[26]. VI. FINAL REMARKS
The Second Peloponnesian War, which began in 431 BCE between the Leagues led by thetwo rival ancient Greek cities of Athens and Sparta, marked for the former – a city which,prior to the war, was living perhaps its most remarkable period (the "age of Pericles") offortunes and great achievements in the arts, culture, and philosophy – a regrettable departurefrom its very recent glorious past. As if the war wasn’t enough a tribulation, Fate bestowedon the city of Athens and neighboring city-states, around 430 BCE, the great plague whichclaimed the lives of about a quarter of its population, including its "first citizen" – the greatPericles (495-429 BCE). Idiosyncratic of the time, the people believed the plague was apunishment from the gods. Legend has it that a delegation was, therefore, dispatched to theoracle of Apollo at Delos to seek guidance on how to stave off the plague, to which they wereinstructed to double the size of the altar of Apollo and make sacrifice on it. The problem,however, was that the altar in question was a cube and the instruction was understood – inline with the traditions of ancient Greek geometry – by all as exactly duplicating the volumeof the altar by means of only two instruments: compass and unmarked straightedge. Thepeople went on to diligently double each side of the altar, only to find that the volume of thenew cubical altar was not twofold but eightfold the original. Many more attempts followedbut the problem stubbornly persisted, but thankfully, the plague eventually worn out withtime.Stated in modern terms, the problem was to solve the equation x = 2 with compass63nd unmarked straightedge. This problem, which came to be known as the Delian problem ,bedeviled for centuries the great minds of Antiquity. As Boyer puts it, "Not until some2000 years later was it recognized that the oracle had sardonically proposed an unsolv-able problem"[29]. The proof that x = 2 cannot be solved with compass and unmarkedstraightedge was to be found in analytic geometry, a new mathematical discipline which wasessentially born out of the insights of two Frenchmen in the early half of the seventeenthcentury: René Descartes (1596-1650) and Pierre de Fermat (1601-1665), with the fecundideas of Nicole Oresme (1323?-1382) and François Viète (1540-1603) – both French as well –playing the role of important precursors. Before the advent of this new discipline, geometryhad been primarily only known for about 2000 years in the form of synthetic geometry (thatis a form of geometry which relies only on axioms and construction methods using compassand straightedge) since the days of Thales of Miletus (ca. 624-548 BCE). But, adopting thealgebraic notations introduced by Viéte and improving them himself, Descartes in his LaGéométrie , first published in 1637, went on to systematically study the relationship betweenalgebra and geometry. This allowed him to translate geometrical problems into algebraicequations and infer some important properties of the geometrical figures he was studying.With such a powerful method in his possession, Descartes even went on to solve the gen-eral n − lines Pappus problem whose solution had even escaped "The Great Geometer" ofAntiquity, Appollonius of Perga (ca. 262 - ca. 190 BCE).The history of mathematics runs abundant with episodes like the above, and we couldmake a long list of them. If any lesson is to be learnt from these, it cannot but be thefact that seemingly intractable mathematical problems can be hacked with just the right setof tools in hand. The process, almost invariably, involves translating the original probleminto an equivalent problem of a different field. In this sense, the HP transformation doesfor multispin magnetic resonance what algebra, so elegantly, did for geometry. And theresults presented in this paper, §IV, are a clear indication that the HP transformation hasall the resemblance of a much needed tool in multispin magnetic resonance. The majorityof these results, for example when we considered in §IV B 1 b the energy levels and theirrespective degeneracy in a system of equivalent spins, are hardly imaginable with the usualspin representations. One important drawback of the results presented above, though, isthat since the dimension of the block matrices B n does not scale polynomially with thenumber of spins N , the simulation of the magnetic resonance spectra of multispins based on64he eigendecomposition of these block matrices is also not going to scale polynomially with N , unless one takes advantage of the eventual presence of groups of equivalent spins (withcardinality ≥ ) in the system, as we saw in §IV B 2. In fact, if we ignore the presence ofthese equivalent groups, the scaling is exponential in the number of spins N , but still farbetter than the exponentially scaling one has with the tout court approach. If we consider amultiset A = { N } with Hamiltonian given by ˆ H o , Eq. (145), for example, we see that thecomputational cost of diagonalizing H o according to the tout court approach scales O (cid:0) Np (cid:1) if the computational cost (according to the algorithm we are using) of diagonalizing a matrixof dimension M is O ( M p ) , where p is a positive real number. With the HP transformation,we will have to diagonalize the matrices B , B , . . . , B N/ (without loss of generality, weassume N is even) to determine the eigenvalues and eigenfunctions of ˆ H o . Knowing thatthe dimension of B n in this case is given by the binomial coefficient (cid:0) Nn (cid:1) , Eq. (100), itthen follows that the computational cost of diagonalizing ˆ H o within the HP transformationframework is: O N/ X n =0 (cid:18) Nn (cid:19) p N ≫ ∼ O Np √ p (cid:18) πN (cid:19) ( p − / ! . (185)Therefore, the complexity of diagonalizing ˆ H o with HP transformation scheme is reducedby a factor proportional to N ( p − / with respect to the tout court approach. With normaldiagonalization algorithms p ≈ , so the reduction factor scales linearly with N . This iscertainly an improvement but not enough when the number of spins exceeds about 15 (andif we don’t exploit the presence of groups of equivalent spins). The argument just outlinedgoes on to show how the diagonalization approach to the simulation of magnetic resonancespectra – despite being exact – is a less favorable one when dealing with a large spin systemwhich lacks groups of equivalent spins with cardinality ≥ . Simulation protocols basedon the HP transformation and a density matrix propagation method like the state-spacerestricted (SSR)[4, 5] method are sure to give a better performance.Density matrix propagation methods for the simulation of multispin magnetic resonancespectra are heavily dependent on the ability to easily index the multispin states at variousstages of the algorithm. This is where the HP transformation and the index compressionmap η o become even more vital. The analytical elements these two mathematical toolsbring to a propagation method can also be advantageous and may help simplify some heavycomputations. 65he multispin Hamiltonians considered in this paper are all given by the sum of operatorswhich are proportional to the q = 0 component of some spherical tensors. Given thatthe leading term in the spin Hamiltonian of a significant number of systems studied inmagnetic resonance (especially in high-resolution problems) are also proportional to the q = 0 component of some spherical tensor, the results presented here are very useful inperturbative calculations on such systems.We have also shown how with the Holstein-Primakoff bosons, we may easily translatecertain problems in multispin magnetic resonance into counting problems. This connectionmay be an indication of a profound and no less consequential relationship between magneticresonance and discrete mathematics; a relationship, which, if ultimately fully understood,might help free – to a significant degree – multispin magnetic resonance from the spell ofthe curse that is of dimensionality. On this note, it is only reasonable that we envisage afuture where the full power and glory of discrete mathematics are unequivocally brought tobear in the study of magnetic resonance. Here, computer algebra – now almost ubiquitousin modern computing languages and which is finding its way also into magnetic resonance[30, 31] – is going to play a central role in the simulation of multispin magnetic resonancespectra. And the recurrent use of generating functions like G A , Ω ( q ) and G A ,λ ( q ) in this paperstrongly suggests the pivotal role polynomial algebra is also going to assume in magneticresonance simulations.It is reported that when the Delian problem was put before Plato (ca. 427 - ca. 347),he interpreted it as a command to the people to study geometry. By the same token, thispaper is a call to integrate discrete mathematics into the theory of magnetic resonance.We hope the discussions and results presented here may kindle at once the curiosity of theReader to further investigate the relation between magnetic resonance, second quantizedspin representations like the HP and Schwinger bosons, and discrete mathematics. Thecurrent paper only minimally scratches the surface of the subject. As Feynman once said,"There’s plenty of room at the bottom". 66 CKNOWLEDGMENTS
The author expresses his gratitude to Prof. Henrik Koch and Prof. Antonino Polimenofor stimulating discussions. [1] P. L. Corio, Chemical Reviews , 363 (1960).[2] T. Holstein and H. Primakoff, Phys. Rev. , 1098 (1940).[3] J. A. Gyamfi and V. Barone, J. Phys. A: Mathematical and Theoretical , 105202 (2018).[4] I. Kuprov, N. Wagner-Rundell, and P. J. Hore, Journal of Magnetic Resonance , 241 (2007).[5] I. Kuprov, Journal of Magnetic Resonance , 45 (2008).[6] R. P. Stanley, “Enumerative combinatorics, vol. i, 2nd edition,” (online; last accessed 14-June-2019).[7] G. H. Fuller, Journal of Physical and Chemical Reference Data , 835 (1976).[8] P. W. Atkins and R. S. Friedman, Molecular Quantum Mechanics (OUP Oxford, 2011).[9] J. Sakurai and J. Napolitano,
Modern Quantum Mechanics (Addison-Wesley, 2011).[10] T. Lancaster, S. Blundell, and S. Blundell,
Quantum Field Theory for the Gifted Amateur (OUP Oxford, 2014).[11] D. Kaiser, B. Chen, D. Derbes, Y. Ting, and R. Sohn,
Quantum Field Theory: Lectures ofSidney Coleman (World Scientific, 2018).[12] R. Zare,
Understanding Spatial Aspects in Chemistry and Physics: Angular Momentum , BakerLecture Series (Wiley, 2018).[13] J. A. Gyamfi and V. Barone, “Magnetic resonance, index compression maps and the holstein -primakoff bosons: Towards a polynomially scaling exact diagonalization of isotropic multispinhamiltonians,” (2018), 1803.10461.[14] L. Dickson,
History of the Theory of Numbers, Volume II: Diophantine Analysis , Dover Bookson Mathematics (Dover Publications, 2005).[15] “Observationes analyticae variae de combinationibus,” From
The Euler Archive – A digital library dedicated to the work and life of Leonhard Euler. http://eulerarchive.maa.org/pages/E158.html (Online; last accessed 14-June-2019).[16] L. Euler, “Various analytic observations on combinations,” (2007), 0711.3656.
17] G. Andrews, “Combinatorics: Ancient & modern,” (OUP Oxford, 2013) Chap. 9, pp. 205–229.[18] G. H. Hardy and S. Ramanujan, Proceedings of the London Mathematical Society s2-17 , 75 (1918).[19] H. Rademacher, Proceedings of the London Mathematical Society s2-43 , 241 (1938).[20] G. E. Andrews,
The theory of partitions , Encyclopedia of Mathematics and its Applications(Cambridge University Press, 1976).[21] A. Zee,
Group Theory in a Nutshell for Physicists , In a Nutshell (Princeton University Press,2016).[22] R. Sachs,
The Physics of Time Reversal (University of Chicago Press, 1987).[23] B. W. Roberts, Philosophy of Science , 315 (2017).[24] E. Wigner, Group Theory: And its Application to the Quantum Mechanics of Atomic Spectra (Academic Press, 1959).[25] J. Jeener, “Superoperators in magnetic resonance,” (Academic Press, 1982) pp. 1 – 51.[26] J. Schwinger,
On Angular Momentum , NYO (Series) (U.S. Atomic Energy Commission, 1952).[27] A. Auerbach,
Interacting Electrons and Quantum Magnetism , Graduate Texts in Contempo-rary Physics (Springer-Verlag New York, 1994).[28] J. Schwinger, Transactions of the New York Academy of Sciences , 170 (1977).[29] C. B. Boyer, Scientific American , 40 (1949).[30] C. Anand, A. Bain, and Z. Nie, Journal of Magnetic Resonance , 200 (2007).[31] X. Filip and C. Filip, Journal of Magnetic Resonance , 95 (2010)., 95 (2010).