# Quantum Derivation of the Bloch Equations Excluding Relaxation

Authors: Eric R. Johnston
Journal: Concepts in Magnetic Resonance Part A, Bridging Education and Research
 (2022)

---

## Abstract

The equation of motion of the density matrix of an ensemble of identical spin-1/2 nuclei subject to a rotating-frame radiofrequency field and Zeeman frequency offset is derived from the Schrodinger equation and shown to be equivalent to the magnetization differential equations originally proposed by Bloch (excluding relaxation). The quantum and classical differential equations are then integrated.


---

## Body

## 1. Introduction

It is well known that the magnetization differential equations (excluding relaxation) of Bloch [1] have a quantum counterpart [2] in the equation of motion of the density matrix under an appropriate Hamiltonian operator. A compact derivation of this equation is presented, and the resulting quantum and classical differential equations are then integrated. The aim throughout is to explicitly illustrate the use of quantum-mechanical principles and matrix methods in formulating and solving this problem. Dirac notation is used to more transparently visualize matrix structures and manipulations. The reader may find it useful to refer first to the overview of Section 3 before beginning Section 2.


## 2. Theory and Results

### 2.1. Derivation of the General Equation of Motion of the Density Matrix

The Schrodinger equation in Dirac notation [3] is(1)ddtψt=−iHψt.H is a time-independent Hermitian operator (square matrix) and ψt is a ket (normalized column vector) state function [4]. The corresponding equation for the adjoint bra ψt (the row vector complex conjugate transpose) is(2)ddtψt=ψtiH,making use of the self-adjoint property H†=H of a Hermitian matrix.For an ensemble of identical particles, the outer product of the ket and bra,(3)ψtψt=σt,is a matrix σt referred to as the density matrix or operator. The time dependence of the density matrix is calculated by differentiating equation (3).(4)ddtσt=ddt|ψt〉〈ψt|=ddt|ψt〉〈ψt|+|ψt〉ddt〈ψt|=−iHσt+σtiH=−iHσt−σtH=−iH,σt,obtained using equations (1), (2), and commutator notation.Equation (4) is known as the Liouville–von Neumann equation of motion of the density matrix.


### 2.2. Form of the Density Matrix

The density matrix elements for a spin−1/2 ensemble are constructed from the expectation values of the spin angular momentum operatorsIx, Iy, or Iz. These are given by the trace (sum of diagonal elements) of the product of the relevant spin angular momentum operator Ix, Iy, or Iz and the density matrix σ and expressed as ensemble projections s¯x, s¯y or s¯zof the appropriate angular momentum in half-integer units [4, 5].(5)Ix=TrIxσ=Tr120110σ11σ12σ21σ22=12σ12+σ21=12s¯x,(6)Iy=TrIyσ=Tr120−ii0σ11σ12σ21σ22=i2σ12−σ21=12s¯y,(7)Iz=TrIzσ=Tr12100−1σ11σ12σ21σ22=12σ11−σ22=12s¯z.Additionally, the trace of the density matrixTrσ=σ11+σ22=1. Solving for individual density matrix elements σ is found to be(8)σ=σ11σ12σ21σ22=121+s¯zs¯x−is¯ys¯x+is¯y1−s¯z.The density matrix is Hermitian with diagonal elements proportional classically to longitudinal z-magnetization and off-diagonal elements proportional to (complex) transverse magnetization [2, 5]. It is idempotent σ2=σ as verified by matrix multiplication using s¯x2+s¯y2+s¯z2=1 for a pure state and Trσ2=1. A geometric formulation of the density matrix is also given in the Appendix. It suffices henceforth to omit the diagonal constant and proportionality factor in equation (8) and express the density matrix as(8a)σ=σ11σ12σ21σ22=IzIx−iIyIx+iIy−Iz,with Ix2+Iy2+Iz2=1.


### 2.3. Form of the Hamiltonian Operator

For a rotating-frame radiofrequency fieldω1=γB1rad/s along the x-axis and a (positive) Zeeman frequency offset Δ=ωrf−ω0rad/s, the Hamiltonian operator H is given by(9)H=ω1Ix+ΔIz=12Δω1ω1−Δ.


### 2.4. Evaluation of the Equation of Motion Equation (4)

Using the density matrixσ of equation (8a) and the Hamiltonian operator H of equation (9) in equation (4), we find(10)ddtσ=ω1Iy−iΔIx−iIy+iω1IziΔIx+iIy−iω1Iz−ω1Iy.Then,d/dtIα=TrIαd/dtσα=x,y,z and(11)ddtIx=Tr120110ddtσ=−ΔIy,(12)ddtIy=Tr120−ii0ddtσ=ΔIx−ω1Iz,(13)ddtIz=Tr12100−1ddtσ=ω1Iy.Equations (11)–(13) may be assembled in the matrix form(14)ddtIxIyIz=0−Δ0Δ0−ω10ω10IxIyIz.As the expectation values are proportional to magnetizations, equations (11)–(13) are seen to be equivalent to the coupled differential equations of Bloch [1]. They may also be written in a compact form d/dtIα=−iIα,H using the cyclic commutation relations Iα,Iβ=iIγ of the spin operators [2].


### 2.5. Superoperator Form of the Liouville–von Neumann Equation

A commutator superoperatorΓ can also be constructed from H using direct products [4].(15)Γ=H⊗E−E⊗H=12Δω1ω1−Δ⊗1001−121001⊗Δω1ω1−Δ=120−ω1ω10−ω12Δ0ω1ω10−2Δ−ω10ω1−ω10.If the elements ofσ in equation (8a) are arrayed as a column supervector σ^, equation (4) becomes(16)ddtσ^=−iΓσ^=−i20−ω1ω10−ω12Δ0ω1ω10−2Δ−ω10ω1−ω10IzIx−iIyIx+iIy−Iz.By solvingd/dtIx,d/dtIy and d/dtIz, we obtain equations (11)–(13).


### 2.6. Integration of the Equation of Motion Equation (4)

The integrated solutions of Schrodinger equations (1), (2) are(17)ψt=exp−iHtψ0,(18)ψt=ψ0expiHt.Accordingly(19)σt=ψtψt=exp−iHtψ0ψ0|expiHt=exp−iHtσ0expiHt.Equation (19) provides a means of calculating the time evolution of the density matrix (unitary transformation) from some initial state with the exponential operator exp−iHt and its adjoint.Equation (4) may be recovered by differentiating equation (19), as follows:(20)ddtσt=−iHσt+exp−iHtσ0iHexpiHt=−iHσt+σtiH=−iH,σt,asH and expiHt commute.


### 2.7. Form of the Exponential OperatorR=exp−iHt

AsIx and Iz do not commute the operator exp−iHt must be explicitly calculated.The diagonal matrix of eigenvalues of the Hamiltonian operatorH of equation (9) is(21)D=12Ω00−Ω,Ω=Δ2+ω121/2.The corresponding normalized eigenvector matrixU is found to be(22)U=cosθ/2−sinθ/2sinθ/2cosθ/2=12,with tanθ=ω1/Δ and expressing U as a row vector of kets.AsU is unitary UT=U−1 and H=UDU−1. The operator R=exp−iHt is then given by(23)R=Uexp−iDtU−1=12c−is1c+is2=c11+22−is11−22=c−iscosθ−issinθ−issinθc+iscosθ=12,using half-angle formulas with c,s=cosΩ/2t,sinΩ/2t and expressing R as a row vector of kets. As R is unitary, the adjoint is given by(24)R†=R−1=expiHt=c+iscosθissinθissinθc−iscosθ=12,using Dirac notation to write R† as the corresponding column vector of bras.


### 2.8. Time Evolution of the Density Matrix for Various Initial States

(a)
With the system initially at equilibriumIz=1 and the density matrix of equation (8a) evolves according to equation (19). Using Dirac notation for R and R†, we may write(25)σt=R100−1R†=11−22=cos2θ+csin2θissinθ+sinθcosθ1−c−issinθ+sinθcosθ1−c−cos2θ−csin2θ,using half-angle formulas and now withc,s=cosΩt,sinΩt (the reader should verify equations (23) and (25)). The trace expressions of equations (5)–(7) then provide the integrated solutions (26)Ix=sinθcosθ1−cosΩt,(27)Iy=−sinθsinΩt,(28)Iz=cos2θ+sin2θcosΩt.(b)
With the system initially alongxIx=1 and σt=R0110R†=12+21. Along yIy=1 and σt=R0−ii0R†=−i12−21. These equations lead, respectively, to sets of expectation values given by column 1 or column 2 of equation (42).(c)
ForIx=sinθ and Iz=cosθ, the ensemble is aligned with the effective field Beff=Ω/γ and σ is time-independent (the reader should verify this assertion).


### 2.9. Radiofrequency Field along they-Axis

For a rf field along they-axis, the Hamiltonian operator H becomes(29)H=ω1Iy+ΔIz=12Δ−iω1iω1−Δ.The eigenvalues are again those of equation (21), and U is found to be(30)U=cosθ/2−sinθ/2isinθ/2icosθ/2=12,from which R is calculated from equation (23) using equation (30) and its adjoint to be(31)R=exp−iHt=Uexp−iDtU†=c−iscosθ−ssinθssinθc+iscosθ=12,using half-angle formulas with c,s=cosΩ/2t,sinΩ/2t and writing R as a row vector of kets.WithIz=1, the density matrix evolves using equation (31) and its adjoint R† according to(32)σt=R100−1R†=11−22=cos2θ+csin2θssinθ−isinθcosθ1−cssinθ+isinθcosθ1−c−cos2θ−csin2θ,using half-angle formulas and now with c,s=cosΩt,sinΩt.Equations (5)–(7) give the integrated solutions(33)Ix=sinθsinΩt,(34)Iy=sinθcosθ1−cosΩt,(35)Iz=cos2θ+sin2θcosΩt.


### 2.10. Integration of the Bloch Differential Equations

Equation (14) may be recast as(36)ddtM=−KM,where M=MxMyMzT is a column vector of magnetizations and(37)K=0Δ0−Δ0ω10−ω10.The integrated solution is(38)Mt=exp−KtM0.The diagonal eigenvalue matrix ofK is D=diag0iΩ−iΩΩ=Δ2+ω121/2, and the normalized eigenvector matrix U that diagonalizes K is found to be(39)U=122sinθcosθcosθ0i−i2cosθ−sinθ−sinθ=123,expressing U as a row vector of kets. As U is unitary, the matrix adjoint U†=U−1, the corresponding column vector of bras. Using Dirac notation, diagonalization of K can be represented as(40)U−1KU=12301iΩ2−iΩ3=0iΩ−iΩ=D,where, e.g., i represents the ith column of U, i is its row adjoint, i|i=1, and i|j=0.K is given by(41)K=UDU−1=1230iΩ−iΩ123=iΩ22−33.The matrixA=exp−Kt is then(42)A=Uexp−DtU−1=1231c−is2c+is3=11+c22+33−is22−33=sin2θ+ccos2θ−scosθsinθcosθ1−cscosθc−ssinθsinθcosθ1−cssinθcos2θ+csin2θ,with c,s=cosΩt,sinΩt (the reader should verify equation (42)).K is antisymmetric so thatKT=−K,exp−KtT=expKt and exp−Ktexp−KtT=E. A is therefore orthogonal and i|i=1,i|j=0.Equation (38) becomes(43)Mt=AM0.For the system initially at equilibrium,M0=001T and Mt is given by column 3 of equation (42). These solutions are those of equations (26)–(28) obtained by integration of the density matrix.


## 2.1. Derivation of the General Equation of Motion of the Density Matrix

The Schrodinger equation in Dirac notation [3] is(1)ddtψt=−iHψt.H is a time-independent Hermitian operator (square matrix) and ψt is a ket (normalized column vector) state function [4]. The corresponding equation for the adjoint bra ψt (the row vector complex conjugate transpose) is(2)ddtψt=ψtiH,making use of the self-adjoint property H†=H of a Hermitian matrix.For an ensemble of identical particles, the outer product of the ket and bra,(3)ψtψt=σt,is a matrix σt referred to as the density matrix or operator. The time dependence of the density matrix is calculated by differentiating equation (3).(4)ddtσt=ddt|ψt〉〈ψt|=ddt|ψt〉〈ψt|+|ψt〉ddt〈ψt|=−iHσt+σtiH=−iHσt−σtH=−iH,σt,obtained using equations (1), (2), and commutator notation.Equation (4) is known as the Liouville–von Neumann equation of motion of the density matrix.


## 2.2. Form of the Density Matrix

The density matrix elements for a spin−1/2 ensemble are constructed from the expectation values of the spin angular momentum operatorsIx, Iy, or Iz. These are given by the trace (sum of diagonal elements) of the product of the relevant spin angular momentum operator Ix, Iy, or Iz and the density matrix σ and expressed as ensemble projections s¯x, s¯y or s¯zof the appropriate angular momentum in half-integer units [4, 5].(5)Ix=TrIxσ=Tr120110σ11σ12σ21σ22=12σ12+σ21=12s¯x,(6)Iy=TrIyσ=Tr120−ii0σ11σ12σ21σ22=i2σ12−σ21=12s¯y,(7)Iz=TrIzσ=Tr12100−1σ11σ12σ21σ22=12σ11−σ22=12s¯z.Additionally, the trace of the density matrixTrσ=σ11+σ22=1. Solving for individual density matrix elements σ is found to be(8)σ=σ11σ12σ21σ22=121+s¯zs¯x−is¯ys¯x+is¯y1−s¯z.The density matrix is Hermitian with diagonal elements proportional classically to longitudinal z-magnetization and off-diagonal elements proportional to (complex) transverse magnetization [2, 5]. It is idempotent σ2=σ as verified by matrix multiplication using s¯x2+s¯y2+s¯z2=1 for a pure state and Trσ2=1. A geometric formulation of the density matrix is also given in the Appendix. It suffices henceforth to omit the diagonal constant and proportionality factor in equation (8) and express the density matrix as(8a)σ=σ11σ12σ21σ22=IzIx−iIyIx+iIy−Iz,with Ix2+Iy2+Iz2=1.


## 2.3. Form of the Hamiltonian Operator

For a rotating-frame radiofrequency fieldω1=γB1rad/s along the x-axis and a (positive) Zeeman frequency offset Δ=ωrf−ω0rad/s, the Hamiltonian operator H is given by(9)H=ω1Ix+ΔIz=12Δω1ω1−Δ.


## 2.4. Evaluation of the Equation of Motion Equation (4)

Using the density matrixσ of equation (8a) and the Hamiltonian operator H of equation (9) in equation (4), we find(10)ddtσ=ω1Iy−iΔIx−iIy+iω1IziΔIx+iIy−iω1Iz−ω1Iy.Then,d/dtIα=TrIαd/dtσα=x,y,z and(11)ddtIx=Tr120110ddtσ=−ΔIy,(12)ddtIy=Tr120−ii0ddtσ=ΔIx−ω1Iz,(13)ddtIz=Tr12100−1ddtσ=ω1Iy.Equations (11)–(13) may be assembled in the matrix form(14)ddtIxIyIz=0−Δ0Δ0−ω10ω10IxIyIz.As the expectation values are proportional to magnetizations, equations (11)–(13) are seen to be equivalent to the coupled differential equations of Bloch [1]. They may also be written in a compact form d/dtIα=−iIα,H using the cyclic commutation relations Iα,Iβ=iIγ of the spin operators [2].


## 2.5. Superoperator Form of the Liouville–von Neumann Equation

A commutator superoperatorΓ can also be constructed from H using direct products [4].(15)Γ=H⊗E−E⊗H=12Δω1ω1−Δ⊗1001−121001⊗Δω1ω1−Δ=120−ω1ω10−ω12Δ0ω1ω10−2Δ−ω10ω1−ω10.If the elements ofσ in equation (8a) are arrayed as a column supervector σ^, equation (4) becomes(16)ddtσ^=−iΓσ^=−i20−ω1ω10−ω12Δ0ω1ω10−2Δ−ω10ω1−ω10IzIx−iIyIx+iIy−Iz.By solvingd/dtIx,d/dtIy and d/dtIz, we obtain equations (11)–(13).


## 2.6. Integration of the Equation of Motion Equation (4)

The integrated solutions of Schrodinger equations (1), (2) are(17)ψt=exp−iHtψ0,(18)ψt=ψ0expiHt.Accordingly(19)σt=ψtψt=exp−iHtψ0ψ0|expiHt=exp−iHtσ0expiHt.Equation (19) provides a means of calculating the time evolution of the density matrix (unitary transformation) from some initial state with the exponential operator exp−iHt and its adjoint.Equation (4) may be recovered by differentiating equation (19), as follows:(20)ddtσt=−iHσt+exp−iHtσ0iHexpiHt=−iHσt+σtiH=−iH,σt,asH and expiHt commute.


## 2.7. Form of the Exponential OperatorR=exp−iHt

AsIx and Iz do not commute the operator exp−iHt must be explicitly calculated.The diagonal matrix of eigenvalues of the Hamiltonian operatorH of equation (9) is(21)D=12Ω00−Ω,Ω=Δ2+ω121/2.The corresponding normalized eigenvector matrixU is found to be(22)U=cosθ/2−sinθ/2sinθ/2cosθ/2=12,with tanθ=ω1/Δ and expressing U as a row vector of kets.AsU is unitary UT=U−1 and H=UDU−1. The operator R=exp−iHt is then given by(23)R=Uexp−iDtU−1=12c−is1c+is2=c11+22−is11−22=c−iscosθ−issinθ−issinθc+iscosθ=12,using half-angle formulas with c,s=cosΩ/2t,sinΩ/2t and expressing R as a row vector of kets. As R is unitary, the adjoint is given by(24)R†=R−1=expiHt=c+iscosθissinθissinθc−iscosθ=12,using Dirac notation to write R† as the corresponding column vector of bras.


## 2.8. Time Evolution of the Density Matrix for Various Initial States

(a)
With the system initially at equilibriumIz=1 and the density matrix of equation (8a) evolves according to equation (19). Using Dirac notation for R and R†, we may write(25)σt=R100−1R†=11−22=cos2θ+csin2θissinθ+sinθcosθ1−c−issinθ+sinθcosθ1−c−cos2θ−csin2θ,using half-angle formulas and now withc,s=cosΩt,sinΩt (the reader should verify equations (23) and (25)). The trace expressions of equations (5)–(7) then provide the integrated solutions (26)Ix=sinθcosθ1−cosΩt,(27)Iy=−sinθsinΩt,(28)Iz=cos2θ+sin2θcosΩt.(b)
With the system initially alongxIx=1 and σt=R0110R†=12+21. Along yIy=1 and σt=R0−ii0R†=−i12−21. These equations lead, respectively, to sets of expectation values given by column 1 or column 2 of equation (42).(c)
ForIx=sinθ and Iz=cosθ, the ensemble is aligned with the effective field Beff=Ω/γ and σ is time-independent (the reader should verify this assertion).


## 2.9. Radiofrequency Field along they-Axis

For a rf field along they-axis, the Hamiltonian operator H becomes(29)H=ω1Iy+ΔIz=12Δ−iω1iω1−Δ.The eigenvalues are again those of equation (21), and U is found to be(30)U=cosθ/2−sinθ/2isinθ/2icosθ/2=12,from which R is calculated from equation (23) using equation (30) and its adjoint to be(31)R=exp−iHt=Uexp−iDtU†=c−iscosθ−ssinθssinθc+iscosθ=12,using half-angle formulas with c,s=cosΩ/2t,sinΩ/2t and writing R as a row vector of kets.WithIz=1, the density matrix evolves using equation (31) and its adjoint R† according to(32)σt=R100−1R†=11−22=cos2θ+csin2θssinθ−isinθcosθ1−cssinθ+isinθcosθ1−c−cos2θ−csin2θ,using half-angle formulas and now with c,s=cosΩt,sinΩt.Equations (5)–(7) give the integrated solutions(33)Ix=sinθsinΩt,(34)Iy=sinθcosθ1−cosΩt,(35)Iz=cos2θ+sin2θcosΩt.


## 2.10. Integration of the Bloch Differential Equations

Equation (14) may be recast as(36)ddtM=−KM,where M=MxMyMzT is a column vector of magnetizations and(37)K=0Δ0−Δ0ω10−ω10.The integrated solution is(38)Mt=exp−KtM0.The diagonal eigenvalue matrix ofK is D=diag0iΩ−iΩΩ=Δ2+ω121/2, and the normalized eigenvector matrix U that diagonalizes K is found to be(39)U=122sinθcosθcosθ0i−i2cosθ−sinθ−sinθ=123,expressing U as a row vector of kets. As U is unitary, the matrix adjoint U†=U−1, the corresponding column vector of bras. Using Dirac notation, diagonalization of K can be represented as(40)U−1KU=12301iΩ2−iΩ3=0iΩ−iΩ=D,where, e.g., i represents the ith column of U, i is its row adjoint, i|i=1, and i|j=0.K is given by(41)K=UDU−1=1230iΩ−iΩ123=iΩ22−33.The matrixA=exp−Kt is then(42)A=Uexp−DtU−1=1231c−is2c+is3=11+c22+33−is22−33=sin2θ+ccos2θ−scosθsinθcosθ1−cscosθc−ssinθsinθcosθ1−cssinθcos2θ+csin2θ,with c,s=cosΩt,sinΩt (the reader should verify equation (42)).K is antisymmetric so thatKT=−K,exp−KtT=expKt and exp−Ktexp−KtT=E. A is therefore orthogonal and i|i=1,i|j=0.Equation (38) becomes(43)Mt=AM0.For the system initially at equilibrium,M0=001T and Mt is given by column 3 of equation (42). These solutions are those of equations (26)–(28) obtained by integration of the density matrix.


## 3. Discussion

The results of Section 2 are summarized here to delineate the steps leading from the Schrodinger eq. (1) to the differential eq. (14) of the expectation values of the spin angular momentum operators and their subsequent integration.(1)
The general differential equation of motion of the density matrix equation (4) is first derived from the Schrodinger equation (1) and its adjoint equation (2).(2)
The density matrix equation (8) for an ensemble of identical spin−1/2 nuclei is then constructed from the expectation values of the spin angular momentum operators.(3)
The Hamiltonian operator equation (9) for a rotating-frame rf field and Zeeman frequency offset is formulated.(4)
The time dependence of the density matrix elements is then calculated (equation (10)). Suitable combinations of these give coupled differential equations (11)–(13) for the expectation values of the spin operators.(5)
An equivalent superoperator formulation of the equation of motion is presented.(6)
Integration of the Schrodinger equations (1), (2) leads to equation (19) describing the time evolution of the density matrix via unitary transformation using an exponential operator and its adjoint.(7)
A diagonalization method is used to calculate the necessary exponential operator (equation (23)).(8)
The time evolution of the density matrix equation (19) is then calculated for the Hamiltonian operator of equation (9) with the initial condition Iz=1, and suitable combinations of density matrix elements give equations (26)–(28) for the spin operator expectation values. Other initial conditions are considered.(9)
A shifted rf field along they-axis (equation (29)) leads to expectation values given by equations (33)–(35) for the system initially at equilibrium.(10)
Finally, a diagonalization method is used to integrate the Bloch differential equation (36), giving equation (43).


## 4. Conclusion

As the expectation values of the operators representing the projection of the spin angular momentum alongx, y or z are proportional to the respective magnetizations, the quantum and classical differential equations of motion (and their integrated solutions) for a spin−1/2 ensemble are shown to be equivalent. It is noteworthy that the Bloch equations, originally proposed phenomenologically [1] using a classical argument, are quantum-mechanical in origin. The equivalence arises from the correspondence [2] between (a) the (classical) cross product coupling of the nuclear magnetic moment vector with the effective field vector and (b) the (quantum-mechanical) commutation relations governing the spin angular momentum operators used to construct the density matrix and the effective field Hamiltonian operator.



---

*Source: 4461529-2022-10-10.xml*