# Solution of the Bloch Equations including Relaxation

Authors: Eric R. Johnston
Journal: Concepts in Magnetic Resonance Part A
 (2020)

---

## Abstract

The magnetization differential equations of Bloch are integrated using a matrix diagonalization method. The solution describes several limiting cases and leads to compact expressions of wide validity for a spin ensemble initially at equilibrium.


---

## Body

## 1. Introduction

In 1949 Torrey used Laplace transforms to provide [1] the first solution of the differential equations proposed by Bloch [2] for the magnetization components of a spin ensemble. His results are somewhat cumbersome and contain some errors. Although the problem is fundamental, a general solution including relaxation does not appear in any of the standard NMR texts, with one partial exception [3]. The problem has been revisited several times employing third-order differential equations [4, 5] and Laplace transforms [6] to give unwieldy solutions using somewhat opaque derivations. The first-order differential equations are directly integrated here using a matrix diagonalization method.


## 2. Bloch Equations and Their Integration

The Bloch equations for a collection of identical spins−1/2 in the frame rotating at ωrf are(1)ddtMxMyMz=−KMxMyMz+00R1M0,where(2)K=R2Δ0−ΔR2ω10−ω1R1.R1andR2 are the longitudinal and transverse relaxation rates in s−1, and Δ=ω0−ωrf and ω1=γB1 are the resonance frequency offset and the rf amplitude for a B1 field along the x-axis, in radians/s. γ is the (positive) gyromagnetic ratio, and M0 is the equilibrium magnetization.Defining a magnetization vectorV=MxMyMzT, the integrated solution of equation (1) is(3)Vt=exp−KtV0−Vss+Vss=Xexp−DtX−1V0−Vss+Vss,where K=XDX−1. The problem is, thus, deducing the roots of K=D, the matrix X and inverse X−1 which diagonalize K and the steady-state magnetization vector Vss.


### 2.1. Evaluation of Roots

K−λ1=0 gives the following cubic equation:(4)λ3−R1+2R2λ2+R22+2R1R2+Δ2+ω12λ−R1R22+R1Δ2+R2ω12=0.Choosing roots of the formλ1=a, λ2,3=b±iΩ gives the corresponding cubic equation:(5)λ−aλ2−2bλ+b2+Ω2=λ3−a+2bλ2+b2+2ab+Ω2λ−ab2+Ω2=0.Equating the last term in equations (4) and (5), and with R2<ω1and/orΔ in (4) and b<Ω in (5),(6)a=R1cos2θ+R2sin2θ≡Rθ,using tanθ=ω1/Δ and Ω≡Δ2+ω121/2.b is obtained by equating the coefficients of the second term in equations (4) and (5):(7)b=12R1+2R2−Rθ=R2−12R2−R1sin2θ.These expressions agree with Torrey ([1], equation 59) and with Abragam ([3]. p. 70).ForR2=1s−1 and ω1and/orΔ=2π×10Hz, for example, the above approach is valid and avoids the explicit solution of the cubic equation (4).


### 2.2. Calculation ofX

X is obtained by evaluating the three cofactors of K−λ1 for λ1=Rθandλ2,3=b±iΩ. Choosing the third row of K−λ1, the cofactors arecofactor 1 =Δω1cofactor 2 =−ω1R2−λcofactor 3 =R2−λ2+Δ2Omitting all (small) relaxation rate difference terms (which are exactly zero forR1=R2) and dividing all elements by ω1 gives X:(8)X=ΔΔΔ0iΩ−iΩΔ2ω1−ω1−ω1.


### 2.3. Calculation ofX−1

X−1 is formed by constructing the matrix of all cofactors of X, taking the transpose, and dividing by the determinant X [7]. The result is(9)X−1=12Ω22ω12Δ02ω1Δ−iΩ−ω1ΔiΩ−ω1.These may be rewritten in a more compact form usingtanθ=ω1/Δ and Ω:(10)X=cosθcosθcosθ0i−icosθtanθ−sinθ−sinθ,(11)X−1=122sinθtanθ02sinθcosθ−i−sinθcosθi−sinθ.Finally, we calculate the matrixA=Xexp−DtX−1:(12)A=Xexp−Rθtexp−btc−isexp−btc+isX−1,wherec=cosΩts=sinΩt.The elements ofA are(13)A11=sin2θexp−Rθt+cos2θexp−btcosΩt,A22=exp−btcosΩt,A12=−A21=−cosθexp−btsinΩt,A23=−A32=−sinθexp−btsinΩt,A13=A31=sinθcosθexp−Rθt−exp−btcosΩt,A33=cos2θexp−Rθt+sin2θexp−btcosΩt.


### 2.4. Steady States

The steady-state magnetizations are found by setting equation (1) to zero and using Cramer’s rule [7]:(14)Mxss=R1Δω1M0d,(15)Myss=−R1R2ω1M0d,(16)Mzss=R1R22+Δ2M0d,where d=R1R22+R1Δ2+R2ω12.ForR2<ω1and/orΔ, these may be simplified using equation (6):(17)Mxss⟶R1RθsinθcosθM0,(18)Myss⟶−R1R2sinθRθΩM0≈0,(19)Mzss⟶R1Rθcos2θM0.The integrated solution equation (3) for the initial condition Mz0=M0 is(20)Mxt=A13M0−A11Mxss−A12Myss−A13Mzss+Mxss,(21)Myt=A23M0−A21Mxss−A22Myss−A23Mzss+Myss,(22)Mzt=A33M0−A31Mxss−A32Myss−A33Mzss+Mzss.In the absence of relaxation, these are the well-known Bloch equations [8, 9]:(23)Mxt=A13M0=sinθcosθM01−cosΩt,(24)Myt=A23M0=−sinθM0sinΩt,(25)Mzt=A33M0=cos2θM0+sin2θM0cosΩt.


## 2.1. Evaluation of Roots

K−λ1=0 gives the following cubic equation:(4)λ3−R1+2R2λ2+R22+2R1R2+Δ2+ω12λ−R1R22+R1Δ2+R2ω12=0.Choosing roots of the formλ1=a, λ2,3=b±iΩ gives the corresponding cubic equation:(5)λ−aλ2−2bλ+b2+Ω2=λ3−a+2bλ2+b2+2ab+Ω2λ−ab2+Ω2=0.Equating the last term in equations (4) and (5), and with R2<ω1and/orΔ in (4) and b<Ω in (5),(6)a=R1cos2θ+R2sin2θ≡Rθ,using tanθ=ω1/Δ and Ω≡Δ2+ω121/2.b is obtained by equating the coefficients of the second term in equations (4) and (5):(7)b=12R1+2R2−Rθ=R2−12R2−R1sin2θ.These expressions agree with Torrey ([1], equation 59) and with Abragam ([3]. p. 70).ForR2=1s−1 and ω1and/orΔ=2π×10Hz, for example, the above approach is valid and avoids the explicit solution of the cubic equation (4).


## 2.2. Calculation ofX

X is obtained by evaluating the three cofactors of K−λ1 for λ1=Rθandλ2,3=b±iΩ. Choosing the third row of K−λ1, the cofactors arecofactor 1 =Δω1cofactor 2 =−ω1R2−λcofactor 3 =R2−λ2+Δ2Omitting all (small) relaxation rate difference terms (which are exactly zero forR1=R2) and dividing all elements by ω1 gives X:(8)X=ΔΔΔ0iΩ−iΩΔ2ω1−ω1−ω1.


## 2.3. Calculation ofX−1

X−1 is formed by constructing the matrix of all cofactors of X, taking the transpose, and dividing by the determinant X [7]. The result is(9)X−1=12Ω22ω12Δ02ω1Δ−iΩ−ω1ΔiΩ−ω1.These may be rewritten in a more compact form usingtanθ=ω1/Δ and Ω:(10)X=cosθcosθcosθ0i−icosθtanθ−sinθ−sinθ,(11)X−1=122sinθtanθ02sinθcosθ−i−sinθcosθi−sinθ.Finally, we calculate the matrixA=Xexp−DtX−1:(12)A=Xexp−Rθtexp−btc−isexp−btc+isX−1,wherec=cosΩts=sinΩt.The elements ofA are(13)A11=sin2θexp−Rθt+cos2θexp−btcosΩt,A22=exp−btcosΩt,A12=−A21=−cosθexp−btsinΩt,A23=−A32=−sinθexp−btsinΩt,A13=A31=sinθcosθexp−Rθt−exp−btcosΩt,A33=cos2θexp−Rθt+sin2θexp−btcosΩt.


## 2.4. Steady States

The steady-state magnetizations are found by setting equation (1) to zero and using Cramer’s rule [7]:(14)Mxss=R1Δω1M0d,(15)Myss=−R1R2ω1M0d,(16)Mzss=R1R22+Δ2M0d,where d=R1R22+R1Δ2+R2ω12.ForR2<ω1and/orΔ, these may be simplified using equation (6):(17)Mxss⟶R1RθsinθcosθM0,(18)Myss⟶−R1R2sinθRθΩM0≈0,(19)Mzss⟶R1Rθcos2θM0.The integrated solution equation (3) for the initial condition Mz0=M0 is(20)Mxt=A13M0−A11Mxss−A12Myss−A13Mzss+Mxss,(21)Myt=A23M0−A21Mxss−A22Myss−A23Mzss+Myss,(22)Mzt=A33M0−A31Mxss−A32Myss−A33Mzss+Mzss.In the absence of relaxation, these are the well-known Bloch equations [8, 9]:(23)Mxt=A13M0=sinθcosθM01−cosΩt,(24)Myt=A23M0=−sinθM0sinΩt,(25)Mzt=A33M0=cos2θM0+sin2θM0cosΩt.


## 3. Results

Limiting forms of equation (3) are discussed in this section.


### 3.1. Case 1: Resonant NutationΔ=cosθ=0,b=R1+R2/2,andMz0=M0

(26)A11=exp−R2t,A12=A21=A13=A31=0,A22=A33=exp−btcosω1t,A23=−A32=−exp−btsinω1t,(27)Mxss=Mxt=0,Myss=−R1ω1M0,Mzss=R1R2ω12M0,(28)Myt=A23M0+Myss1−A22−A23Mzss,(29)Mzt=A33M0+Mzss1−A33−A32Myss.My and Mz interconvert at rate ω1 and decay to (small) steady states.


### 3.2. Case 2: Free Precession/Mz Relaxation ω1=sinθ=0andMx0=Mzss=M0

(30)A11=A22=exp−R2tcosΔt,A12=−A21=−exp−R2tsinΔt,A13=A31=A23=A32=0,A33=exp−R1t,(31)Mxt=A11M0=M0exp−R2tcosΔt,Myt=A21M0=M0exp−R2tsinΔt,Mzt=−A33M0+Mzss=M01−exp−R1t.Mx and My interconvert at rate Δ and decay to zero as Mz returns to equilibrium.


### 3.3. Case 3: Spin-Locked RelaxationMx0=sinθM0,Mz0=cosθM0,andMyss≈0

(32)Mxt=A11sinθM0+A13cosθM0−A11Mxss−A13Mzss+Mxss,Myt=A21sinθM0+A23cosθM0−A21Mxss−A23Mzss,Mzt=A31sinθM0+A33cosθM0−A31Mxss−A33Mzss+Mzss.Employing the steady states of equations (17)–(19), these become(33)Mxt=sinθM0exp−Rθt+R1RθsinθcosθM01−exp−Rθt,Myt=0,Mzt=cosθM0exp−Rθt+R1Rθcos2θM01−exp−Rθt.The magnetization vector relaxes to a steady state along the effective fieldBeff=Ω/γ [10, 11].


### 3.4. Case 4: General SolutionMz0=M0

Equations (20)–(22) may be recast using the steady states of equation (17)–(19) as(34)Mxt=A13M0+Mxss1−exp−Rθt,Myt=A23M0,Mzt=A33M0+Mzss1−exp−Rθt.These expressions agree with equations (60)–(62) of Torrey [1] and with Abragam [3].For a weak rf fieldω1<Δcosθ⟶1,sin2θ⟶0, they reduce to(35)Mxt=ω1ΔM01−exp−R2tcosΔt,Myt=−ω1ΔM0exp−R2tsinΔt,Mzt≈M0,in agreement with Slichter ([12], p. 35).


### 3.5. Case 5: Equal Relaxation RatesRθ=b=R,Mz0=M0

An exact solution of the Bloch equations is given by equations (20)–(22) using the full steady-state expressions (14)–(16). Using the steady states of equation (17)–(19), they become(36)Mxt=sinθcosθM01−exp−RtcosΩt,Myt=−sinθM0exp−RtsinΩt,Mzt=cos2θM0+sin2θM0exp−RtcosΩt.


## 3.1. Case 1: Resonant NutationΔ=cosθ=0,b=R1+R2/2,andMz0=M0

(26)A11=exp−R2t,A12=A21=A13=A31=0,A22=A33=exp−btcosω1t,A23=−A32=−exp−btsinω1t,(27)Mxss=Mxt=0,Myss=−R1ω1M0,Mzss=R1R2ω12M0,(28)Myt=A23M0+Myss1−A22−A23Mzss,(29)Mzt=A33M0+Mzss1−A33−A32Myss.My and Mz interconvert at rate ω1 and decay to (small) steady states.


## 3.2. Case 2: Free Precession/Mz Relaxation ω1=sinθ=0andMx0=Mzss=M0

(30)A11=A22=exp−R2tcosΔt,A12=−A21=−exp−R2tsinΔt,A13=A31=A23=A32=0,A33=exp−R1t,(31)Mxt=A11M0=M0exp−R2tcosΔt,Myt=A21M0=M0exp−R2tsinΔt,Mzt=−A33M0+Mzss=M01−exp−R1t.Mx and My interconvert at rate Δ and decay to zero as Mz returns to equilibrium.


## 3.3. Case 3: Spin-Locked RelaxationMx0=sinθM0,Mz0=cosθM0,andMyss≈0

(32)Mxt=A11sinθM0+A13cosθM0−A11Mxss−A13Mzss+Mxss,Myt=A21sinθM0+A23cosθM0−A21Mxss−A23Mzss,Mzt=A31sinθM0+A33cosθM0−A31Mxss−A33Mzss+Mzss.Employing the steady states of equations (17)–(19), these become(33)Mxt=sinθM0exp−Rθt+R1RθsinθcosθM01−exp−Rθt,Myt=0,Mzt=cosθM0exp−Rθt+R1Rθcos2θM01−exp−Rθt.The magnetization vector relaxes to a steady state along the effective fieldBeff=Ω/γ [10, 11].


## 3.4. Case 4: General SolutionMz0=M0

Equations (20)–(22) may be recast using the steady states of equation (17)–(19) as(34)Mxt=A13M0+Mxss1−exp−Rθt,Myt=A23M0,Mzt=A33M0+Mzss1−exp−Rθt.These expressions agree with equations (60)–(62) of Torrey [1] and with Abragam [3].For a weak rf fieldω1<Δcosθ⟶1,sin2θ⟶0, they reduce to(35)Mxt=ω1ΔM01−exp−R2tcosΔt,Myt=−ω1ΔM0exp−R2tsinΔt,Mzt≈M0,in agreement with Slichter ([12], p. 35).


## 3.5. Case 5: Equal Relaxation RatesRθ=b=R,Mz0=M0

An exact solution of the Bloch equations is given by equations (20)–(22) using the full steady-state expressions (14)–(16). Using the steady states of equation (17)–(19), they become(36)Mxt=sinθcosθM01−exp−RtcosΩt,Myt=−sinθM0exp−RtsinΩt,Mzt=cos2θM0+sin2θM0exp−RtcosΩt.


## 4. Discussion

Equation (3) describes a number of experimental situations.


### 4.1. Case 1: Resonant Nutation

Resonant nutation (equations (28) and (29) ) was described by Torrey in his original paper [1]. The effective relaxation rate is the average of the longitudinal and transverse relaxation rates [3].


### 4.2. Case 2: Free Precession and Relaxation

In the absence of an rf field, the transverse components interconvert and relax to zero (the free induction decay) as the longitudinal component, initially zero, relaxes independently to equilibrium (equation (31)).


### 4.3. Case 3: Spin-Locked Relaxation

Orientation of the magnetization vector parallel to the effective field suppresses precession and results in a single-exponential approach to equilibrium, affording the longitudinal and transverse relaxation rates using equations (6) and (33) [10, 11].


### 4.4. Case 4: General Solution

Equation (34) presents in compact form the solutions originally given by Torrey [1] and by Morris and Chilvers [6] as Laplace expressions and the tabulations of Madhu and Kumar [4, 5] for a spin ensemble initially at equilibrium. They are valid provided(37)R2s−1<Ωrads−1,which holds for most cases of practical interest. In the example given by Madhu and Kumar [4, 5],(38)R1=10s−1,R2=100s−1,Δ=ω1=2π×1kHz.Accordingly, from equations (6), (7), and (17)–(19),(39)Rθ=55s−1,b=77.5s−1,Ω=2π×1414Hz,Mxss=Mzss=M011,Myss=−M0110×2π≈0.These values do not appear to reproduce the figures presented in [5].


### 4.5. Case 5: Equal Relaxation Rates

The solutions (20)–(22) for equal relaxation rates are exact provided the full steady states of equation (14)–(16) are used. The inequality of Case 4 leads to the simpler expressions (36). We note also that setting ω1=2Δ results in Rθ=b=R1+2R2/3.


### 4.6. Neglect of Relaxation

For rf amplitude and precession terms which are large compared to relaxation rates equations (23)–(25) pertain. They are useful, for example, in describing selective (on-resonance) excitation with (off-resonance) signal suppression [9] as in the following example (using Hz units).


#### 4.6.1. On-Resonance Rotation

(40)cosθ⟶0,v1tp=14,Mytp=−M0,Mxtp=Mztp=0.


#### 4.6.2. Off-Resonance Rotation

(41)ΩtP=1,MztP=M0,MxtP=MytP=0,Ωtp=1requiresΔtp=154sotp=154Δsec.Then,v1=1/4tp=Δ/15Hz.For a 5 KHz offsettp=193.6μsec and v1=1291Hz.The on-resonance magnetization is rotated to the –y axis by the rf pulse, whereas the off-resonance magnetization undergoes an excursion that returns it to the z-axis.


## 4.1. Case 1: Resonant Nutation

Resonant nutation (equations (28) and (29) ) was described by Torrey in his original paper [1]. The effective relaxation rate is the average of the longitudinal and transverse relaxation rates [3].


## 4.2. Case 2: Free Precession and Relaxation

In the absence of an rf field, the transverse components interconvert and relax to zero (the free induction decay) as the longitudinal component, initially zero, relaxes independently to equilibrium (equation (31)).


## 4.3. Case 3: Spin-Locked Relaxation

Orientation of the magnetization vector parallel to the effective field suppresses precession and results in a single-exponential approach to equilibrium, affording the longitudinal and transverse relaxation rates using equations (6) and (33) [10, 11].


## 4.4. Case 4: General Solution

Equation (34) presents in compact form the solutions originally given by Torrey [1] and by Morris and Chilvers [6] as Laplace expressions and the tabulations of Madhu and Kumar [4, 5] for a spin ensemble initially at equilibrium. They are valid provided(37)R2s−1<Ωrads−1,which holds for most cases of practical interest. In the example given by Madhu and Kumar [4, 5],(38)R1=10s−1,R2=100s−1,Δ=ω1=2π×1kHz.Accordingly, from equations (6), (7), and (17)–(19),(39)Rθ=55s−1,b=77.5s−1,Ω=2π×1414Hz,Mxss=Mzss=M011,Myss=−M0110×2π≈0.These values do not appear to reproduce the figures presented in [5].


## 4.5. Case 5: Equal Relaxation Rates

The solutions (20)–(22) for equal relaxation rates are exact provided the full steady states of equation (14)–(16) are used. The inequality of Case 4 leads to the simpler expressions (36). We note also that setting ω1=2Δ results in Rθ=b=R1+2R2/3.


## 4.6. Neglect of Relaxation

For rf amplitude and precession terms which are large compared to relaxation rates equations (23)–(25) pertain. They are useful, for example, in describing selective (on-resonance) excitation with (off-resonance) signal suppression [9] as in the following example (using Hz units).


### 4.6.1. On-Resonance Rotation

(40)cosθ⟶0,v1tp=14,Mytp=−M0,Mxtp=Mztp=0.


### 4.6.2. Off-Resonance Rotation

(41)ΩtP=1,MztP=M0,MxtP=MytP=0,Ωtp=1requiresΔtp=154sotp=154Δsec.Then,v1=1/4tp=Δ/15Hz.For a 5 KHz offsettp=193.6μsec and v1=1291Hz.The on-resonance magnetization is rotated to the –y axis by the rf pulse, whereas the off-resonance magnetization undergoes an excursion that returns it to the z-axis.


## 4.6.1. On-Resonance Rotation

(40)cosθ⟶0,v1tp=14,Mytp=−M0,Mxtp=Mztp=0.


## 4.6.2. Off-Resonance Rotation

(41)ΩtP=1,MztP=M0,MxtP=MytP=0,Ωtp=1requiresΔtp=154sotp=154Δsec.Then,v1=1/4tp=Δ/15Hz.For a 5 KHz offsettp=193.6μsec and v1=1291Hz.The on-resonance magnetization is rotated to the –y axis by the rf pulse, whereas the off-resonance magnetization undergoes an excursion that returns it to the z-axis.


## 5. Conclusion

The differential equations (1) of Bloch [2] are integrated with a matrix diagonalization method to give the solution equation (3). It correctly describes a number of experimental situations including resonant nutation, free precession and relaxation, and spin-locked relaxation. Equation (3) is exact for the case of equal longitudinal and transverse relaxation rates and leads to the general equation (34) for a spin ensemble initially at equilibrium.



---

*Source: 8819956-2020-09-02.xml*