Magnetic Resonance, Index Compression Maps and the Holstein-Primakoff Bosons: Towards a Polynomially Scaling Exact Diagonalization of Isotropic Multispin Hamiltonians
aa r X i v : . [ c ond - m a t . o t h e r] O c t IOP/123-QED
Magnetic Resonance, Index Compression Maps and the Holstein-Primakoff Bosons:Towards a Polynomially Scaling Exact Diagonalization of Isotropic MultispinHamiltonians.
J.A. Gyamfi a) and V. Barone Scuola Normale Superiore di Pisa, Piazza dei Cavalieri 7, 56126 Pisa,Italy. (Dated: October 17, 2018)
Matrix diagonalization has long been a setback in the numerical simulation of themagnetic resonance spectra of multispin systems since the dimension of the Hilbertspace of such systems grows exponentially with the number of spins – a problemcommonly referred to as the "curse of dimensionality". In this paper, we proposetwo mathematical instruments which, when harmoniously combined, could greatlyhelp surmount to a fair degree and in a systematic manner the curse of dimension-ality. These are: 1) the Holstein-Primakoff bosons and 2) what we have termed the"index compression maps". These two allow a bijective mapping of (multi)spin statesto integers. Their combination leads to the block diagonalization of the multispinHamiltonian, thus a computationally exact way of diagonalizing the latter but whichalso reduces significantly the computational cost. We also show that the eigenvectorsand eigenvalues of the Liouvillian operator can be easily determined once those of therelated multispin Hamiltonian are known. Interestingly, the method also enables ananalytical characterization of the multispin Hilbert space – a feat hardly attainablewith other approaches.We illustrate the method here by showing how a general static isotropic multi-spin Hamiltonian could be exactly diagonalized with very less computational cost.Nonetheless, we emphasize that the method could be applied to study numerousquantum systems defined on finite Hilbert spaces and embodied with at most pair-wise interactions. a) Email corresponding author: jerryman.gyamfi@sns.it . INTRODUCTION It is of common knowledge that the computational cost of simulating many-body quan-tum systems on classical computers grow exponentially with the dimension D of the system’sHilbert space, making the computation hardly feasible when D exceeds some threshold .Such computational headaches fall under what has been christened the "curse of dimension-ality".When it comes to simulating the magnetic resonance spectrum (say EPR or NMR) of achemical species, the problem is particularly felt and somehow frustrating because, on topof the exponential growth of the computational cost, one ends up with even more sparsematrices as the number of nonzero spins in the system increases.The relentless growth of the sparseness of the matrices is an indication that the dynamicsof these systems are effectively encoded in certain subspaces of the Hilbert space, whosedimensions in most cases scale polynomially with the number of particles. Thus, overcomingthe curse of dimensionality is, in the long run, equivalent to finding an efficient way ofdetermining the said subspaces and restricting to them the numerical simulation of thesystem’s dynamics. In magnetic resonance, computational protocols like the adaptive state-space restriction (SSR) method proposed by Kuprov et. al. , for example, constructs aneffective subspace by excluding entanglements beyond a certain fixed order, thus yieldingan effective subspace which scales polynomially in the total number of spins N . We showhere that it is possible to devise a computationally exact diagonalization scheme for staticisotropic multispin Hamiltonians which reduces significantly the computational cost, withoutexcluding any states or interactions. The method hints at the possibility of a polynomiallyscaling numerical protocol for the simulation of multispin dynamics, and yet more analyticalin nature.We are of the view that the curse of dimensionality can be fairly dealt with if one couldeliminate superfluous mathematical operations during numerical simulations. But this re-quires adequate analytical tools. We may say the mathematical operation at the core of thereason why the computational cost of multispin numerical simulations grows exponentiallywith the number of spins is the Kronecker product. It would thus be most desirable if wehave a set of analytical instruments which could help us "navigate" confidently through thismathematical maze. It should therefore not be of any surprise if central to the method we2resent here are: 1) the Holstein-Primakoff bosons, and 2) what we have termed the "rowindex compression map" η . Together, they expedite the determination of new good quan-tum numbers – besides allowing an easy monitoring of the various spin transitions in thesystem. In addition, η releases us from the burden of carrying along long strings of indexesduring theoretical derivations, making them manageable – not to mention its usefulness innumerical simulations. The great benefit of the index compression maps we discuss belowis that they enable us to construct the Kronecker product at the local level, which presentsgreat advantages. This will be clear very soon.Though we focus here on multispin systems and magnetic resonance, the method couldbe easily employed in other problems, provided the Hilbert space is finite and presents atmost pairwise interactions.The paper is organized as follows: After a brief review of the notations we have adopted(§II), we first consider the index compression map theorem and its consequences in §III. Thisis followed by a brief review of the Holstein-Primakoff transformation in §IV. With thesetwo instruments in hand, we illustrate how they could be used to efficiently diagonalize ina computationally exact manner static isotropic multispin Hamiltonians with less computa-tional cost compared to traditional methods in §V. In §VI, we apply our results to brieflydiscuss the diagonalization of the isotropic static multispin Hamiltonian of the deuteratedhydroxymethyl radical ( • CH OD).
II. NOTATIONS AND TERMINOLOGY
Consider a multiset A of N spins: A := { j , j , . . . , j N } ≡ { j N α α } , where α ∈ { , , . . . , σ } .The index α certainly runs over the distinct elements (in total σ ) of A , and N α is themultiplicity of j α . Notably, N = P σα =1 N α . By distinct elements we mean SU (2) spins withdifferent maximal weights – other characteristics of the spins (like charge, magnetic momentetc.) are not of any merit whatsoever here. For example, the multiset A = {
12 9 , } indicatesany aggregate of nine spin- / and three spin- ; the actual composition could consist of, forinstance, 1) nine electrons and three protons, or 2) five N , four muons – both spin- / s –and three N (spin- ) , etc. In the following, we reserve the Latin letter i to lower index A ’s elements when the latter is understood as A = { j , j , . . . , j N } , while we use the Greekindex α when we mean distinct elements of A .3perators 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 .Finally, we shall refer to nonzero submatrices as simply "submatrices". III. THE INDEX COMPRESSION MAPS, η AND ξ . The Hilbert space H of a given multispin system A is constructed by taking the directproduct of the Hilbert space H i of each individual spins. The basis (in the uncoupledrepresentation) of H is determined in similar fashion. The index compression maps (inparticular, the row index one) we explore below allow one to determine the composition ofeach basis element of the global Hilbert space in terms of the basis elements of the separateHilbert spaces H i . In addition, we gain a bijective mapping between H ’s basis and the setof natural numbers. This is particularly useful as we shall see in the next section.We present below the index compression map theorem in general terms. Its applicationand usefulness go far beyond the applications we propose in this paper. Theorem 1 (The index compression map theorem) . Given the finite set of finite matrices { A , A , . . . , A N } – where the dimension of the matrix A i is m i × n i , ( m i , n i ∈ N ) – and theordered Kronecker product C = A ⊗ · · · ⊗ A N − ⊗ A N (1a) c ν,µ = a ν ,µ · · · a ν N − ,µ N − a ν N ,µ N , (1b) then, there is a one-to-one correspondence between the index ν ( µ ) and the indexes { ν , . . . , ν N − , ν N } ( { µ , . . . , µ N − , µ N } ) given by the mapping η ( ξ ), such that η : G r ⊂ N × N × · · · × N | {z } N −→ N ν = η ( ν , . . . , ν N ) = 1 + N X i =1 W R,i ( ν i − (2)4 : G c ⊂ N × N × · · · × N | {z } N −→ N µ = ξ ( µ , . . . , µ N ) = 1 + N X i =1 W C,i ( µ i − (3) where, W R,i ≡ δ N,i + (1 − δ N,i ) N − − i Y l =0 m N − l (4a) W C,i ≡ δ N,i + (1 − δ N,i ) N − − i Y l =0 n N − l . (4b)The entries of the matrix C resulting from the Kronecker product given in Eq. (1a) areproducts which have as factors an entry from each factor matrix on the RHS. The mapping ξ , termed the "column index compression map", assigns a unique column integer to the(ordered) collection of column indexes of the factors. Similarly, η maps a string of orderedrow indexes of the factors to a unique row index of C . G r ( G c ) is the set of all possibleordered combinations of row (column) indexes of the factor matrices. It can be proved thatboth η and ξ are bijective, thus with inverses. That is, given ν , for example, one can findthe corresponding string of ordered row indexes of the factor matrices.Take for example the following Kronecker product: C = A ⊗ B ⊗ D (5a) c ν,µ = a ν ,µ b ν ,µ d ν ,µ (5b)where, A = a , a , a , a , a , a , a , a , a , B = b , b , b , D = d , d , d , d , . (6)We have here desisted from denoting the factor matrices as A , A and A in order to simplifythe notations, but it is clear that the matrices A, B, D in Eq. (5a) correspond to A , A and A in Eq. (1a), respectively. The map η relates any given set of row indexes ν , ν , ν of the matrices A, B, D , respectively, with a specific row index ν of C . The map ξ does an5nalogous thing but with the column indexes µ , µ , µ and µ . Now, from Eqs. (2) - (4), itfollows that for the particular Kronecker product given in Eq. (5a), the maps η ( ν , ν , ν ) and ξ ( µ , µ , µ ) are respectively of the form, ν = η ( ν , ν , ν ) = 1 + 6( ν −
1) + 2( ν −
1) + ( ν − , (7a) µ = ξ ( ν , ν , ν ) = 1 + 2( µ −
1) + 2( µ −
1) + ( µ − . (7b)Say we want to know which element of C corresponds to the product a , b , d , . In otherwords, we want to know ν and µ given that ν = 2 , µ = 1 (since the element of A con-tributing to the product is a , ), ν = 3 , µ = 1 (given that we have b , ), and ν = 2 , µ = 2 .Inserting these integers into Eqs. (7), we find that ν = η (2 , ,
2) = 12 and µ = ξ (1 , ,
2) = 2 .Hence, a , b , d , = c , .In the other way round, we may have ν and µ of c ν,µ and may want to know whichcomponent of A, B, D gives rise to that specific element of C . For example, we may askwhat are the elements of A, B, D whose product gives rise to the element c , of C . Thiscorresponds to solving the Diophantine equations η ( ν , ν , ν ) = 1 + 6( ν −
1) + 2( ν −
1) + ( ν − , (8a) ξ ( µ , µ , µ ) = 1 + 2( µ −
1) + 2( µ −
1) + ( µ − . (8b)for ν , ν , ν , µ , µ , µ . The dimensions of the factor matrices A, B, D serve as constraintson the unknown variables. A quick strategy in solving these Diophantine equations is thefollowing: Consider for example Eq. (8a). Find the integer among the possible values of ν which maximizes the sum of the first two terms; the sum must not be greater than ν , in thiscase . Find then the integer in the range of ν which maximizes the first three terms of Eq.(8a), with ν equal to the integer found in the previous step. Here too, the sum must notexceed the value of ν . Continue the process until you get to ν N , which must be such that itsvalue, together with the those of ν , . . . , ν N − found earlier, gives exactly ν according to Eq.(2). The same strategy applies in the case of the column index µ . As for Eq. (8a), we findthat ν = 2 , ν = 2 and ν = 1 . Analogously, we find that for Eq. (8b) to hold, the followingmust be true: µ = 3 , µ = µ = 1 . Thus, c , = a , b , d , , which can be easily verified. Forthe procedure outlined to give the right solution, it is important that the summands of η and ξ are written in the same order as done on the RHS of Eqs. (8a)-(8b), i.e. the constantterm must appear first, then the ν (or µ ) term, and so on. To understand why this is so,6t is useful to realize that the compression maps η and ξ may be interpreted as functionswhich convert string of numbers from a "multi-numeral" base to the decimal base (for now,the integer zero is not an allowed digit).It is important to note that the order of the factor matrices determine the equation forthe maps η and ξ . Indeed, a permutation of the factor matrices effects a permutation ofthe elements of the resulting matrix. For example, if instead of the product in Eq. (5a) weconsider the following product C ′ = B ⊗ A ⊗ D (9a) c ′ ν ′ ,µ ′ = b ν ′ ,µ ′ a ν ′ ,µ ′ d ν ′ ,µ ′ (9b)where with respect to Eq. (5a) we have swapped the matrices B and A , we note that themappings η and ξ transform as ν ′ = η ( ν ′ , ν ′ , ν ′ ) = 1 + 6( ν ′ −
1) + 2( ν ′ −
1) + ( ν ′ − , (10a) µ ′ = ξ ( ν ′ , ν ′ , ν ′ ) = 1 + 6( µ ′ −
1) + 2( µ ′ −
1) + ( µ ′ − . (10b)With Eq. (9a), the product a , b , d , = c ′ , . Nonetheless, a , b , d , is still c ′ , .As one can see, the lowest value of both ν and µ as defined in Eqs. (2) and (3) is . Thesame applies to their corresponding strings of factor matrix indexes ν i and µ i . We emphasizethat one can set arbitrarily this lowest value to the integer z , whence η and ξ become η z and ξ z , respectively, and so η z , ξ z : Z × Z × · · · × Z | {z } N −→ Z (11a) ν = η z ( ν , . . . , ν N ) = z + N X i =1 W R,i ( ν i − z ) (11b) µ = ξ z ( µ , . . . , µ N ) = z + N X i =1 W C,i ( µ i − z ) . (11c)In this article, we shall adopt the following compression maps: η and ξ , i.e. we set thelowest value of the indices ν, ν i , µ, µ i to zero.We may apply these results to the construction of finite tensor spaces. Say we havea collection of finite vector spaces {V , . . . , V N } from which we create the vector space V through the ordered tensor product V ≡ V ⊗ · · · ⊗ V N − ⊗ V N . (12)7he basis of V is given by the ordered Kronecker product of the basis of the spaces {V , . . . , V N } . Assume we number progressively (starting from ) the elements of thebasis of each space V i , and then represent each basis element by a column vector. The basis B i of V i may thus be represented, symbolically, by a column vector, B i = e ( i )1 ... e ( i ) d i (13)where e ( i ) k is the assigned k − th basis element of the vector space V i , and d i = dim V i . Then,the following corollary follows from the index compression map theorem: Corollary 1.1.
Let B be the column representation of the basis of V , Eq. (12) , then B = B ⊗ · · · ⊗ B N − ⊗ B N (14a) e ν = e (1) ν ⊗ · · · ⊗ e ( N − ν N − ⊗ e ( N ) ν N (14b) where e ν is the ν − th element of the basis of V , and ν = η ( ν , . . . , ν N ) = 1 + N X i =1 W R,i ( ν i − . (15)Note that in writing Eqs. (13) and (14) we have ignored the column indexes since they arealways equal to 1. At this point, it is crucial we draw the reader’s attention to the following:from Eq. (14) we see that to every tensor product of basis, we can associate a uniqueordered string of the type { ν , . . . , ν N − , ν N } , which in turn corresponds to the unique basis e ν , where ν is given by (15). There is thus a one-to-one correspondence between the string { ν , . . . , ν N − , ν N } and ν , as already stated by the index compression map theorem (theorem1). What this means is that we can forgo working directly with the tensor product of basisand work instead with ν ∈ { , . . . , W R, } and its corresponding string { ν , . . . , ν N − , ν N } .The index ν in Eq. (15) is said to be a "reducible" index, while the indexes { ν , . . . , ν N − , ν N } are its "irreducible" indexes (or components). It is worth noting that W R, , Eq. (4a),coincides with the dimension of the tensor space V . Note however that if, instead of the η ,we employ the map η in Eq. (15), then ν ∈ { , , . . . , W R, − } and every irreducible indexshall have as minimum instead of . 8he application of the index compression mappings in finite tensor Hilbert spaces isstraightforward. We recall here again that in the following we shall be using the map η .The reason why will be explained shortly. IV. THE HOLSTEIN-PRIMAKOFF TRANSFORMATION
As it is well-known, in the uncoupled representation, the basis of the i − th spin Hilbertspace H i is {| j i , m i i} , where j i is the spin quantum number and m i an eigenvalue of theoperator ˆ J i,z . Certainly, the dimension of H i is given by the number of possible values m i can take, which is j i + 1 . Since the basis {| j i , m i i} is a countable set, we may choose toassign an integer to each element of the set. There are numerous possible ways to do sobut let us settle on the following: we start with the basis element with the highest value of m i , i.e. | j i , j i i and assign to it the integer . In other words, we make the transformation | j i , j i i 7→ | i . We then take the next basis element | j i , j i − i and assign to it the integer . So | j i , j i − i 7→ | i , and so forth – up to | j i , − j i i 7→ | j i i . In general, we have effectedthe transformation | j i , m i i 7→ | j i − m i i . The new basis {| i , . . . , | j i i} is the integer repre-sentation of H i ’s basis and the integers define the irreducible indexes on H i . What we havejust accomplished is nothing less than a second quantization. Precisely, it is part of what istermed in the literature as the "Holstein-Primakoff transformation".In concrete terms, the Holstein-Primakoff (HP) transformation is a one-flavored map-ping of spin- j operators to boson creation and annihilation operators. By "one-flavored"we mean only one type of boson is involved. Consider a collection A of N spins: A = { j , j . . . , j N } . Then, the HP transformation for the i − th spin operators reads ˆ J − i ˆ b † i (cid:16) j i ˆ i − ˆ b † i ˆ b i (cid:17) / (16a) ˆ J + i (cid:16) j i ˆ i − ˆ b † i ˆ b i (cid:17) / ˆ b i (16b) ˆ J zi j i ˆ i − ˆ b † i ˆ b i (16c)where as usual, ˆ b i | ν i i = √ ν i | ν i − i , ˆ b † i | ν i i = √ ν i + 1 | ν i + 1 i , ˆ b † i ˆ b i | ν i i = ν i | ν i i , h ˆ b i , ˆ b † i ′ i = δ i,i ′ ˆ i (17)and ν i ∈ { , , , . . . , j i } . The bosons in discussion here are the so-called Holstein-Primakoffbosons (or magnons), and one may think of the spins as vertexes of a complete graph – where9ach vertex can accommodate not more than a certain fixed number of these particles .For this reason, we shall sometimes refer to the spins as "spin vertexes" when working inthe HP spin representation. Moreover, noteworthy is the observation that the occupationnumber of the HP bosons uniquely specifies the orientation of the spin along the axis ofquantization. That is, independent of whether one is dealing with half-integral spins or not,any spin orientation is indicated by a nonnegative integer in the HP transformation. Thisis no trivial computational advantage.The transformations listed in Eq. (16) ensure that the usual spin commutation relationsare obeyed, namely h ˆ J αi , ˆ J βi ′ i = iǫ αβγ δ i,i ′ ˆ J γi , (18)where ǫ αβγ is the three dimensional Levi-Civita symbol.There is thus a one-to-one correspondence between the basis {| j i , m i i} and the integerbasis {| ν i i} under the HP transformation . In explicit terms, | j i , m i i 7→ p ( j i − m i )! (cid:16) ˆ b † i (cid:17) j i − m i | i = | ν i i (19)where ν i = j i − m i .For what concerns magnetic resonance, the N spins may consist of N el unpaired electronsand N nu nuclei of a given radical, for example. Then the (spin) Hilbert space H of the systemmay be defined according to the ordered Kronecker product H = H ⊗ · · · ⊗ H N − ⊗ H N , (20)where H i represents the i − th spin’s Hilbert space. Certainly, a viable basis B for H consistsof the Kronecker product (maintaining the same order as in Eq. (20)) of the basis of theisolated H i : B = {| j , m i ⊗ | j , m i ⊗ · · · ⊗ | j N , m N i} , (21)where − j i ≤ m i ≤ j i . This is obviously the uncoupled representation. Naturally, if we moveto the integer representation, we may associate in a bijective manner a unique integer (re-ducible index, ν ) to each basis element of B , Eq. (21), employing the row index compressionmap η (Corollary 1.1) . Thus, we have | ν i = | ν i ⊗ · · · ⊗ | ν N − i ⊗ | ν N i ≡ | ν , . . . , ν N − , ν N i , (22)10here | ν i i is the integer representation of that basis element of H i which contributes to thetensor product defining | ν i . Having chosen the mapping η , it follows from Eq. (11) that ν = η ( ν , . . . , ν N ) = N X i =1 W R,i ν i , (23)and ν i ∈ { , , . . . , j i } . It should now be clear why the choice η : with it, the enumerationof the normal spin basis elements coincides exactly with the possible occupation numbers ofthe HP bosons.It must be emphasized that while the irreducible indexes { ν i } indicate occupation num-bers of the HP bosons, the reducible ones { ν } enumerate the various possible permutationsof the occupation numbers of these bosons. V. THE ROW INDEX COMPRESSION MAP, THE HP BOSONS ANDTHE DIAGONALIZATION OF STATIC ISOTROPIC MULTISPINHAMILTONIANS
We show in this section, as a matter of illustration, how the mapping η together withthe HP bosons can lead to the block diagonalization of an isotropic static multispin Hamil-tonian ( i.e. one with only isotropic coupling constants, g − factors included), and thus to adiagonalization scheme which is computationally more affordable. In the following, all basesmentioned are to be understood as orthonormal.We now subject the spin multiset A to a magnetic field of intensity B applied along the z − axis, which is chosen as the quantization axis. The isotropic spin Hamiltonian ˆ H spin ofthe system in the uncoupled representation is then ˆ H spin = B X i µ i g i ˆ J zi + X i>i ′ T i,i ′ ˆ J i · ˆ J i ′ (24)where ˆ J i (cid:16) = ˆ J xi + ˆ J yi + ˆ J zi (cid:17) , µ i and g i are the spin operator, the appropriate magneton(Bohr or nucleus) and the (isotropic) g − factor of the i − th spin, respectively, and T i,i ′ isthe coupling constant (for example, the dipole-dipole coupling constant plus others) of theinteraction between spin i and i ′ . ˆ H spin may be written explicitly in matrix representationas an operator G spin acting on the identity operator I ⊗ N of H , i.e. H spin = G spin I ⊗ N (25)11ith, I ⊗ N := N O i =1 I i = I i =1 ⊗ · · · ⊗ I i = N (26a) G spin := B X i G (cid:16) I i , µ i g i ˆ J zi (cid:17) + X i,i ′ ( >i ) X β ∈{ x,y,z } T i,i ′ G (cid:16) I i , ˆ J βi ; I i ′ , ˆ J β ′ i ′ (cid:17) , (26b)where ˆ J βi is the spin angular momentum operator of the i − th spin along the axis β , andwhere G is the substitution operator defined operationally as G ( x ′ , x ; . . . ; x ′ n , x n ) f ( x , . . . , x n ) = f ( x ′ , . . . , x ′ n ) . (27)The steady increase in ˆ H spin ’s sparseness is evident in Eq. (25), if we notice that G spin comprises substitution operators which substitute at most two identity matrices of I ⊗ N withother matrices, independent of the total number of spins N . A. Block diagonalizing the spin Hamiltonian
For starters, we rewrite ˆ H spin according to the HP transformation: ˆ H spin = B X i µ i g i (cid:16) j i − ˆ b † i ˆ b i (cid:17) + X i>i ′ T i,i ′ (cid:16) j i − ˆ b † i ˆ b i (cid:17) (cid:16) j i ′ − ˆ b † i ′ ˆ b i ′ (cid:17) + 12 X i>i ′ T i,i ′ (cid:20) ˆ b † i (cid:16) j i − ˆ b † i ˆ b i (cid:17) / (cid:16) j i ′ − ˆ b † i ′ ˆ b i ′ (cid:17) / ˆ b i ′ + h.c. (cid:21) . (28)Consider now the matrix element D λ (cid:12)(cid:12)(cid:12) ˆ H spin (cid:12)(cid:12)(cid:12) ν E between two arbitrary reducible indexes, λ and ν , of H . It follows from Eq. (28) that D λ (cid:12)(cid:12)(cid:12) ˆ H spin (cid:12)(cid:12)(cid:12) ν E = δ λ ,ν " B µ B N X i =1 g i ( j i − ν i ) + X i>i ′ T i,i ′ ( j i − ν i )( j i ′ − ν i ′ ) + 14 X i>k T i,k ( Y l = i,k δ λ l ,ν l ! δ λ i ,ν i ± δ λ k ,ν k ± h − sgn ( λ i − ν i ) sgn ( λ k − ν k ) i × s ± ( ν i ) s ± ( ν k ) ) (29)where s ± ( ν i ) := s(cid:18) ν i + 1 ∓ (cid:19) (cid:18) j i − ν i + 1 ± (cid:19) . (30)12otice that the matrix element h λ | H spin | ν i reduces to only one of the two terms in Eq.(29): it reduces to the first term if the two reducible indexes λ and ν coincide, otherwise ityields the second term.Regarding the second term, one observes that the necessary and sufficient condition(assuming T i,k = 0 ) for any of its summand to be non nonzero is if the following conditionsare both met: • if of the N spin vertexes, there is no variation in the number of bosons occupying N − spin vertexes ( i.e. we have N − invariant irreducible indexes) as the transitionfrom the initial reducible index ν to the final index λ is made; • if of the remaining two spin vertexes, the transition from ν to λ is such that one losesa boson to the other.These observations actually suggest a good quantum number in the chosen representation;namely, the sum of the irreducible components of the reducible indexes constitute good quan-tum numbers . In other words, the total number of HP bosons are conserved during transi-tions. Indeed, the total boson occupation number operator ˆ N = P Ni =1 ˆ b † i ˆ b i commutes with ˆ H spin . Therefore, a necessary but not sufficient condition for D λ (cid:12)(cid:12)(cid:12) ˆ H spin (cid:12)(cid:12)(cid:12) ν E to be nonzero isthat n ( ν ) = n ( λ ) (31a) n ( ν ) := D ν (cid:12)(cid:12)(cid:12) ˆ N (cid:12)(cid:12)(cid:12) ν E = N X i =1 ν i . (31b)Say N the set { n ( ν ) } , i.e. the set of all possible eigenvalues of the operator ˆ N . Since n is a good quantum number of ˆ N , it follows that all reducible indexes with the same n forma subspace of the total spin Hilbert space H . And any two subspaces characterized by twodifferent n values are orthogonal to each other. It is evident that N ≡ { , , . . . , J } , (32)where J = P i j i . Consequence of these considerations is the following proposition: Proposition 1.
Given a collection of N spins with spin quantum numbers j , . . . , j N whoseinteractions are described by an isotropic static spin Hamiltonian of the form in Eq. (29) , he block-diagonalized H spin consists of exactly J + 1 submatrices, i.e. H spin = J M n =0 B n = diag ( B , B , . . . , B J ) . (33)The complexity of finding the eigenvalues and eigenvectors of H spin thus reduces to thatof finding the same for J + 1 smaller matrices. In particular, if U spin and D spin are H spin ’seigenvector and eigenvalue matrices, respectively, then U spin = J M n =0 U n and D spin = J M n =0 D n (34)where U n and D n are the matrices composed of the eigenvectors and eigenvalues of the n − thsubmatrix in Eq. (33), respectively.It is worth noting from Eq. (33) that the total number of submatrices B n is determinedby the total spin J and not by the total number of spins N . This translates into theimplication that any arbitrary collection of spins having the same total spin number J described by the Hamiltonian ˆ H spin of Eq. (29) will have the same number of submatriceswhen the latter is block-diagonalized, though their dimension may differ from system tosystem. B. Dimension of the submatrices. Density and sparseness of H spin The use of HP bosons allow an easy analytical characterization of the total spin Hilbertspace H . Since all reducible indexes with the same number of HP bosons n span thesame subspace of H , the dimension of the submatrix B n , Ω A ,n , is simply the number ofcompositions of the integer n according to Eq. (31b). We recall that the irreducible indexes ν i in the mentioned equation are subject to the constraint ≤ ν i ≤ j i .The analytical determination of Ω A ,n is an enumerative combinatoric problem . In par-ticular, it can be proved that the generating function, G A , Ω ( q ) , for the integers { Ω A ,n } is ofthe form : G A , Ω ( q ) = σ Y α =1 (cid:16) [2 j α + 1] q (cid:17) N α = J X n =0 Ω A ,n q n , (35)where [ m ] q simply indicates the q − analogue of the integer m , viz. [ m ] q := 1 − q m − q = 1 + q + q + . . . + q m − . (36)14he dimension of the submatrices may be computed straightforward analytically by resortingto the following formula Ω A ,n = X P σα =1 (2 j α +1) s α ≤ n ≤ s α ≤ N α ( − s + ... + s σ (cid:18) N + n − − P σα =1 (2 j α + 1) s α N − (cid:19)(cid:18) N s (cid:19) . . . (cid:18) N σ s σ (cid:19) . (37)We remark that the coefficients { Ω A ,n } are reciprocal and unimodal, namely, Ω A ,n = Ω A , J − n , ∀ n ∈ N (38a) Ω A ,n ≤ Ω A ,n ′ , for n < n ′ and ≤ { n, n ′ } ≤ ⌊ J ⌋ , (38b)respectively. In addition, max { Ω A ,n } = Ω A , ⌊ J ⌋ , where ⌊•⌋ is the floor function.Interestingly, we note from Eq. (37) that Ω A ,n – being the sum of products of binomialcoefficients – implies that the computational cost of diagonalizing the { B n } submatricesis still exponential in N but scales lesser with the latter as compared to diagonalizing thewhole Hamiltonian at once. As an example, consider spin − / s described by a static ˆ H spin of the form in Eq. (29). The dimension of its global Hilbert space is = 1024 ,which means finding its eigenvalues and eigenvectors would require diagonalizing a matrixof dimension . But using the row index compression map η and the HP bosons, thesame eigenvectors and eigenvalues can be determined by diagonalizing eleven (in reality onlynine) matrices, the largest being of dimension . In particular, submatrices of dimension , , , and will appear in pairs, and only one will be of dimension . If we nowadd another spin- / to the system, we will have to diagonalize now a matrix of dimension but employing the scheme explained above we will need to diagonalize smaller matricesof dimension not greater than .Still on the topic of analytic characterization of the multispin Hamiltonian, we may studyfor example the density of H spin . If we define the density ζ ( A ) of the matrix A as the ratiobetween its number of nonzero elements and the total number of entries, then it can beproved that given a collection of N spins and assuming a multispin Hamiltonian of the form(24), the density of H spin has an upper bound given by the expression sup ζ ∈Z [ ℘ ( H spin )] ζ ( H spin ) = P J n =0 Ω A ,n Q σα =1 (2 j α + 1) N α , (39)where ℘ ( H spin ) is the order of H spin and Z [ ℘ ( H spin )] is the set of all possible ζ s of H spin .In the limit case whereby the spin multiset A consists of only spin-1/2s, i.e. A = n N o , one15erives from Eq. (39) that ζ ( H spin ) ≤ − N (cid:18) NN (cid:19) , (40)which for very large N reads approximately ζ ( H spin ) . √ πN . (41)The density of the spin Liouvillian L spin on the other hand has an upper bound of sup ζ ∈Z [ ℘ ( L spin )] ζ ( L spin ) = (cid:16)P J n =0 Ω A ,n (cid:17) Q σα =1 (2 j α + 1) N α . (42)which is the square of that of H spin .The parameter ζ is a good measure of the fraction of the total Hilbert (or Liouville)space which effectively encodes the whole spin dynamics. A closely related measure is thesparseness χ := 1 − ζ . C. Exploiting T -symmetry Although the spin Hamiltonian, Eq. (24), we are concerned with here is real, it fails tobe time-reversal (or T − symmetry) invariant (in the conventional sense) due to the presenceof the Zeeman term. Indeed, as it is well-known, the presence of an external magnetic fieldin a non-spinless system always leads to broken T − symmetry . Nonetheless, the subspaces B n and B J − n are related through T − symmetry. We show below how this relation couldbe exploited to further reduce the computational cost of the numerical simulation of themagnetic resonance spectra of multispin systems.Let ˆ T i represent the time-reversal operator acting on the i − th spin. Then, in the integerrepresentation, the irreducible index | ν i i is transformed by ˆ T i as follows ˆ T i | ν i i = | j i − ν i i ≡ | ν i i , (43)which means the unitary part of ˆ T i (defined up to a phase factor) rotates the {| ν i i} vectorsby an angle of π around the x − axis. The vector | ν i i is said to be the hole complement of | ν i i . As a matter of fact, in the HP representation, the time-reversal symmetry could beinterpreted as the particle-hole symmetry.Since ˆ T i | ν i i = | ν i i (44)16e have that ˆ T i = ˆ i = ⇒ ˆ T i = ˆ T † i .For a multispin system, the total time-reversal operator reads ˆ T = Y i ˆ T i (45)and ˆ T | ν i = ˆ T | ν , ν , . . . , ν N i = | j − ν , j − ν , . . . , j N − ν N i = | W R, − − ν i := | ν i . (46) | ν i is clearly the hole complement of | ν i . Naturally, ˆ T = ˆ T † .We now investigate how ˆ T rotates ˆ H spin . To begin with, it can be easily proved that ˆ T i ˆ b † i ˆ b i ˆ T † i = 2 j i − ˆ b † i ˆ b i , (47a) ˆ T i ˆ J ± i ˆ T † i = ˆ J ∓ i . (47b)Then, as a consequence of the transformations in Eqs. (47), we are led to the followingimportant identity, ˆ H spin − ˆ T ˆ H spin ˆ T † = 2 B X i µ i g i (cid:16) j i − ˆ b † i ˆ b i (cid:17) . (48)That is, the spin Hamiltonian ˆ H spin and its time-reversally rotated counterpart differ bytwice the Zeeman term, which is perfectly in line with what we should expect since thespin-spin interaction terms in ˆ H spin are T − invariant.To understand the importance of Eq. (48), we must bear in mind that if | λ i and | ν i areelements of the basis which span the subspace B n , then their complements (cid:12)(cid:12) λ (cid:11) and | ν i alsospan B J − n . It follows from Eq. (48) that D λ (cid:12)(cid:12)(cid:12) ˆ H spin (cid:12)(cid:12)(cid:12) ν E − D λ (cid:12)(cid:12)(cid:12) ˆ H spin (cid:12)(cid:12)(cid:12) ν E = δ λ,ν B X i µ i g i ( j i − ν i ) , (49)which greatly simplifies the determination of the matrix elements of B J − n if those of B n are known. In concrete terms, B J − n and B n differ only in their diagonal elements, thedifference given by twice the corresponding Zeeman terms of B n . In terms of matrixrepresentation, we have B J − n = T B n T − Z n ∀ n ∈ N (50)17here Z n is the matrix representation of the Zeeman term in Eq. (24) for the subspace B n ,and T in Eq. (50) could be understood, without loss of generality, as the unit antidiagonalmatrix of order Ω A ,n . A straight-away computation of the eigenvalues and eigenvectors of B J − n when those of B n are known could be particularly interesting and useful. D. Eigenvalues and eigenvectors of the Liouvillian ˆˆ L spin It mostly happens that in theoretical derivations, working in Liouville space appearsto be more convenient. We may say that the Liouville space is also ideal for simulatingmagnetic resonance spectra, but workers in the field are well acquainted with the fact thatfrom a practical point of view it is not since the dimension of the Liouville space is that ofthe corresponding Hilbert space squared, which means the computational cost gets roughlysquared with respect to that in the Hilbert space if one tries to diagonalize the Liouvillianfor subsequent calculations. This hurdle can be easily overcome as we explain below.We recall that the matrix representation of the Liouvillian superoperator ˆˆ L spin corre-sponding to ˆ H spin is L spin = H spin ⊗ I H − I H ⊗ H spin , (51)where I H is the matrix representation of H ’s identity operator. It can easily be proved thatthe supermatrix U spin which diagonalizes L spin is related to the eigenvector matrix of H spin , U spin , Eq. (34), as U spin = U spin ⊗ U spin , (52)and the diagonal matrix of L spin , D spin , is related to D spin , Eq. (34), through the identity D spin = D spin ⊗ I H − I H ⊗ D spin . (53)These last two relations are valid in general for real (Hermitian) Hamiltonians, independentof whether the Hamiltonian is time-dependent or not. For complex Hermitian Hamiltonians,they slightly change .Therefore, one understands that if the eigenvalues and eigenvectors of H spin are known,those of the corresponding Liouvillian matrix L spin can be directly computed without havingto diagonalize the superoperator de novo. The significance of the relations (52) and (53) liesin the fact that the computational cost of finding the eigenvalues and eigenvectors of ˆˆ L spin is now essentially the same as finding those of ˆ H spin . So, with a diagonalization protocol18f H spin like the one presented in this article whose computational cost scales very lessexponentially with N , Eqs. (52) and (53) allow simulation in Liouville space without muchextra computational cost. This point is of great consequence because if we consider againthe case of spin- / s, the Liouvillian L spin is a square matrix of dimension , , but through Eqs. (52) and (53), its eigenvectors and eigenvalues can be determined bydiagonalizing a series of matrices ( in total) of dimension , , , , and usingthe scheme presented above. An additional relevance of the Eqs. (52) and (53) has to dowith the fact that they may allow an efficient perturbative computation of the resonancespectrum of the spin system in Liouville space inasmuch as they drastically reduce thecomputational cost of diagonalizing L spin .We conclude this section pointing out that there is a one-to-one correspondence betweenthe total occupation numbers { n } and the set of all possible total spin projections along z (or weights), i.e. { M n } . In particular, M n = J − n . Since the occupation numbersof the HP bosons relate to spin orientations, saying the number of HP bosons is conservedimplies the total spin weights are also conserved: thus ∆ M = 0 , for a unique M n – hence weare dealing with zero-quantum transitions. The observation just made seems to disparagethe importance of the diagonalization scheme presented above in relation to the study ofmultispin dynamics but that should not be the case for the following reason: ˆ H spin as given inEq. (29) may be viewed as the multispin Hamiltonian in some convenient frame of reference,for instance, the molecular frame. Then, upon diagonalizing it, as we have shown above,one can perform a time-dependent (possibly stochastic) rotation of the equations describingthe dynamics to that of the laboratory frame. And this will inevitably generate a dissipativeterm which will induce the transitions observed experimentally. The Hamiltonian in the non-dissipative term of the rotated equation of motion will have the same eigenvalues as ˆ H spin ,but the new eigenvectors will be simply rotated with respect to the original ones and can beeasily determined. The latter may be used as the basis in some perturbative calculation of thecontribution of the dissipative term to the evolution of the system under some appropriateassumptions. The rotation may also be directly performed in Liouville space effortlessly.The important thing noteworthy here is that the procedure we just described begins withthe diagonalization of the static spin Hamiltonian in the chosen convenient frame. We shallreturn to this procedure in a future article.19 I. APPLICATION: THE • CH OD RADICAL, A CASE STUDY.
Suppose we want to simulate the electron paramagnetic resonance (EPR) spectrum of thedeuterated hydroxymethyl radical ( • CH OD). We will have to write its spin Hamiltonianin some basis. We assume the experimental conditions are such that the isotropic spinHamiltonian given in Eq. (24) is suitable, with all relevant constants known. Our aim hereis not to work out the final spectrum but to show the usefulness of the mapping η and theHP bosons.We have four species of nonzero spin in the • CH OD radical: the unpaired electron (spin- / ), the two methyl hydrogen atoms (each of spin- / ) and the deuterium isotope (spin- ).The spin multiset A we are considering here is therefore A = n
12 3 , o . And the Hilbert spaceof the radical is thus H = H ( e ) ⊗ H (1) ⊗ H (2) ⊗ H ( D ) , (54)where H ( e ) is the Hilbert space of the unpaired electron, while H (1) , H (2) and H ( D ) representthe Hilbert space of the two methyl hydrogens and deuterium, respectively. In writing (54)we have chosen a specific order, which is of course not unique. Nonetheless, we must stickto the chosen order when constructing H ’s basis in the uncoupled representation.Now, in the integer representation, we know that all uncoupled states | j i , m i i undergothe transformation | j i , m i i 7→ | j i − m i i (55)where the nonnegative integer ( j i − m i ) is the number of HP bosons occupying that spinvertex. For example, for the spin- / particles we have | / , / i 7→ | i , | / , − / i 7→ | i , (56)while for the spin- species we have | , i 7→ | i , | , i 7→ | i , | , − i 7→ | i . (57)The Hilbert space H of the radical is of dimension , and each of its basis elements isspecified by the tensor product involving one basis element of each of its subspaces. H ’s20asis may be represented conveniently as (following the same order as in Eq. (54)), | ν ) z }| { | | ... | = | ν i z }| { | i| i ⊗ | ν i z }| { | i| i ⊗ | ν i z }| { | i| i ⊗ | ν i z }| { | i| i| i . (58)where for the sake of clarity, we have indicated the basis elements of H in the integerrepresentation as | ν ) . It is clear from Eq. (58) that to each of the kets on the left hand sideis associated a unique tensor product of the basis elements of the subspaces, thus the indices | ν ) of H are reducible. For example, |
1) = | i ⊗ | i ⊗ | i ⊗ | i = | , , , i . In matrix terms,that corresponds to the following column vector: |
1) = ⊗ ⊗ ⊗ (59)As mentioned earlier, the reducible indexes { ν } enumerate the various possible occupationsof the HP bosons on the graph of spin vertexes. For instance, |
1) = | , , , i relates to theHP boson configuration whereby all the spin- / vertexes are empty and only the spin − vertex is occupied by a single HP boson.The row-index compression map η which relates the string { ν , ν , ν , ν } to ν for thechosen order of the Kronecker product in Eq. (54) is ν = η ( ν , ν , ν , ν ) = 12 ν + 6 ν + 3 ν + ν . (60)From this equation, one can determine the reducible index ν for any given string ν , ν , ν , ν ,and vice versa. For example, one can easily determine that |
15) = | i ⊗ | i ⊗ | i ⊗ | i .Having assumed that the Hamiltonian is of the form in Eq. (24), we are immediately ledto conclude that the total occupation number operator ˆ N of the HP bosons commutes withthe spin Hamiltonian. Thus, reducible indexes | ν ) containing the same number of HP bosonsform a subspace orthogonal to all others. Since the total spin of the radical is J = 5 / , theremust be J + 1 = 6 orthogonal subspaces B n , each characterized by the total number of21osons n ∈ { , , , , , } each basis element of the subspace contains. The dimension Ω A ,n of these subspaces is readily determined, for example, by means of the generating function G A , Ω ( q ) , Eq. (35). In the case of • CH OD, we have G A , Ω ( q ) = [2] q · [3] q = X n =0 Ω A ,n q n = 1 + 4 q + 7 q + 7 q + 4 q + q . (61)So, instead of diagonalizing the matrix H spin , which is of dimension , to obtain its eigen-values and eigenvectors, we can alternatively accomplish the same by diagonalizing foursmaller matrices: two of dimension and the other two of dimension . The subspace B and B J are always of dimension one, and correspond to multispin states whereby all theindividual spins are in their maximum weight states, i.e. m i = j i ∀ i (for B ) and lowestweight (for B J ), respectively.The basis of the subspaces { B n } may be determined by writing all the possible strings ofirreducible indexes: i.e. {| , , , i , | , , , i , . . . , | , , , i , | , , , i} , and then group-ing them according to the total number of HP bosons n they contain. As we haveseen above, strings containing the same n HP bosons span the subspace B n . An evenmore convenient alternative in determining the basis is to exploit the mapping η , Eq.(60), and the particle-hole symmetry, Eq. (46). For instance, the basis of B are cer-tainly {| , , , i , | , , , i , | , , , i , | , , , i} ≡ {| , | , | , | } . Then, due to theparticle-hole symmetry, we know that the complements of B ’s basis, i.e. {| , | , | , | } ,span B . The subspaces B n for • CH OD, their dimension and basis elements are reportedin Table I.Once we know the basis elements of the various subspaces, we can construct their matrixrepresentations through Eq. (29). To save computational time, we can make use of Eq.(50), or equivalently (49). The eigenvalues and eigenvectors of H spin then follow from thediagonalization of these smaller matrices. In addition, the density ζ , Eq. (39), of theradical’s spin Hilbert space cannot exceed . . In other words, the percentage of H which effectively encodes the dynamics of the system cannot be greater than ∼ of H .The eigenvectors and eigenvalues of the radical’s Liouvillian operator can be easily computedfrom Eqs. (52) and (53), respectively. 22 ubspace B n , n ( ν ) = P i =1 ν i Dimension of subspace Basis elements, | ν ) B | B | , | , | , | B | , | , | , | , | , | , | B | , | , | , | , | , | , | B | , | , | , | B | Table I:
Orthogonal subspaces B n and their dimensions for the radical • CH OD, assuminga spin Hamiltonian of the generic form given in Eq. (24) . The integer representation | ν ) ofthe basis elements spanning each subspace depends on the order chosen when constructingthe tensor space H , Eq. (54). VII. CONCLUSION
We have shown that through the row index compression map η and the HP transforma-tion, one can easily generate good quantum numbers (in general, integers defined by somecomposition law) which characterize different subspaces of the Hilbert space. The greatmerit of this technique is that neither the matrix H spin nor L spin needs to be constructed inorder to perform the diagonalization: all one needs are the submatrices B n of H spin . Butthese can also be easily constructed on account of Eqs. (29) and (31). The procedure is asfollows: 1) group together all reducible indexes | ν i containing the same total number of HPbosons n , 2) construct the submatrix B n for each group using Eq. (29), and 3) diagonalizethe obtained submatrices.Naturally, the presented block diagonalization scheme of the multispin Hamiltonianthrough the use of η and the HP bosons reduces exponentially the computational complex-ity of the problem of finding the eigenvalues and eigenvectors of the latter. We have alsoshown that the same eigenvectors and eigenvalues can be used to derive the correspondingeigenvectors and eigenvalues of the corresponding Liouvillian superoperator ˆˆ L spin of ˆ H spin ,thus saving us from the even higher computational cost of diagonalizing the former all overagain. Furthermore, we emphasize that this diagonalization scheme should be seen as the23rst step towards a more robust, efficient and analytical procedure to overcome the curseof dimensionality, possibly scaling polynomially with N . We mention that once the subma-trices have been obtained, one can employ other efficient numerical protocols to diagonalizethe larger ones. This will further minimize the computational cost.Besides allowing an easy, yet analytic, monitoring of the various multispin states andtransitions (because η , like all η c , is bijective) – not to mention the yoke of carrying alonglong strings of indexes which it frees us from – another key advantage of the mapping η has to do with the fact that it indubitably allows an analytic characterization of multi-spin Hamiltonians without much effort as we saw in sections V A and V B where we wereable to predict how many subspaces we should be expecting and their dimensions. We callthe reader’s attention to the fact that such results were easily obtained due to the con-nection the mapping η and HP bosons create between spin composition and enumerativecombinatorics .That spin Hamiltonians of the form in Eq. (24) decompose into submatrices is per se notany news. There are a number of ways to prove the same. Corio, for example, has extensivelytreated the topic in his seminal book . What the literature lacked, as far as we know, is ageneral, comprehensive and systematic way to effectively realize such decomposition. Andthat is what the present article aims at.The procedure we have presented can be easily extended to anisotropic multispin Hamil-tonians. There, the inherent symmetries of the coupling tensors together with other sym-metries may be exploited in determining the good quantum numbers in the integer repre-sentation. Moreover, the method presented also lends itself to the study of various quantumsystems defined on finite Hilbert spaces, these include a variety of spin lattice models andmany variants of the simple symmetric exclusion processes (see for example ), to say theleast.Finally, we remark that the presence of square roots of operators in the HP transformationmay prove to be inconvenient in certain scenarios. By resorting to the Schwinger bosons ,one avoids these square roots at the cost of dealing with two flavors. The use of Schwingerbosons is worth pursuing as it might probably lead to an even more efficient diagonalizationscheme. 24 CKNOWLEDGMENT
The authors are very grateful to Prof. Antonino Polimeno (Univ. of Padova) and Prof.Giorgio J. Moro (Univ. of Padova) for their discerning comments upon reading an abridgedversion of the manuscript. JAG would also like to thank Mr. Andrea Piserchia (SNS) for themany interesting discussions they have had over the past months. The authors employedcomputer facilities at SMART@SNS laboratory, to which they are greatly grateful. Theresearch leading to these results has received funding from the European Research Councilunder the European Union’s Seventh Framework Program (FP/2007-2013) / ERC GrantAgreement n. [320951].
REFERENCES Auerbach, A.,
Interacting Electrons and Quantum Magnetism , Graduate Texts in Con-temporary Physics (Springer-Verlag New York, 1994). Barone, V., Gyamfi, J. A., and Piserchia, A., “ Electron paramagnetic resonance: Volume25,” (Royal Society of Chemistry, 2016) Chap. 4, pp. 98–156. Corio, P. L.,
Structure of High-Resolution NMR Spectra (Academic Press, 1966). Feynman, R. P., International Journal of Theoretical Physics , 467 (1982). Fuller, G. H., Journal of Physical and Chemical Reference Data , 835 (1976). Gyamfi, J. A. and Barone, V., J. Phys. A: Mathematical and Theoretical , 105202 (2018). Haake, F.,
Quantum Signatures of Chaos , Physics and astronomy online library (Springer,2001). Holstein, T. and Primakoff, H., Phys. Rev. , 1098 (1940). Kuprov, I., Wagner-Rundell, N., and Hore, P. J., J. Magn. Reson. , 241 (2007). Mendonça, J. R. G., Journal of Physics A: Mathematical and Theoretical , 295001 (2013). In the general case, the Liouvillian L of the complex Hamiltonian H defined on the finiteHilbert space H , is L = H ⊗ I H − I H ⊗ H T , (62)25here H T is the transpose of H . If U ( U ) and D ( D ) are the eigenvector and eigenvaluematrices, respectively, of H ( L ) , then U = U ⊗ (cid:0) U − (cid:1) T , D = D ⊗ I H − I H ⊗ D ..