Consider a collection of water molecules.  Each proton in each molecule is a spin ½ particle with a magnetic moment m = γS, γ = gnqe/mp (SI units).
If we place a sample in a uniform magnetic field B = (0,0,B0), then, neglecting the interaction between different particles, the Hamiltonian of each proton is H = -m∙B = -mzB0.  (Here we treat the space coordinates macroscopically.)
H = -γB0Sz = -ω0Sz.
The eigenstates of H are the common eigenstates of S2 and Sz, {|+>, |->}.  The eigenvalues are E + = -ω0ħ/2,  E- = +ω0ħ/2.

Assume a collection of water molecules, placed in uniform magnetic field B = (0,0,B0), is in thermal equilibrium with its surroundings.  We then are dealing with a statistical mixture of states.  The density operator for the system is ρ = p+|+><+| + p-|-><-|.
p± = N exp(-E±/(kT)) = (1/Z) exp(-E±/(kT)).
We have Tr(ρ) = <+|ρ|+> + <-|ρ|-> = (1/Z)(exp(ω0ħ/(2kT)) + exp(-ω0ħ/(2kT)) = 1.
The partition function Z is therefore given by Z = 2 cosh(ω0ħ/(2kT)).
The density matrix at t = 0 in the {|+>, |->} basis is given by
.
How does the density matrix evolve?
(d/dt)ρ(t) = (1/(iħ))[H,ρ(t)]
The matrix of H in the {|+>, |->} basis is

We therefore have
11/dt = 0,  dρ22/dt = 0,  dρ12/dt = iω0ρ12,  dρ21/dt = -iω0ρ21.
This yields

We can now calculate <Sx>, <Sy>, and <Sz> for our sample.
<Sx> = Tr(ρSx).

Similarly, <Sy> = 0.
<Sz> = Tr(ρSxz).
.
The macroscopic magnetization M of the water sample is given by M = nγ<S>, where n is the number of protons per unit volume.  We therefore have M = (0,0,Mz) where Mz = nγ<Sz>.

Assume that at t = 0 a time-varying field B1(t) is added to B0.  Assume that for t > 0 we have B = (B1cosωt, -B1sinωt, B0).  The z-component of B is constant, while the component of B perpendicular to the z-axis rotates cw about the z-axis with frequency ω.
Assume that B1 << B0.  The proton Hamiltonian is now given by
H(t) = -γSB(t) = -γB0Sz - γB1(cosωt Sx - sinωt Sy) = -ω0Sz - ω1(cosωt Sx - sinωt Sy),
with ω1 = γB1.

Consider the state vector |ψ(t)> describing a single proton, |ψ(t)> = U(t,0)|ψ(0)>.
Let R(φ) describe a ccw rotation about the z-axis and let U(R) be the rotation operator for the ccw rotation R(φ).
U(R) = exp(-(i/ħ)Szφ).
We have
|ψ(t)> = U(t,0)|ψ(0)>,  U(R)|ψ(t)> = U(R)U(t,0)|ψ(0)>.
U(R) can be viewed as an operator that rotates every state vector ccw through an angle φ or an operator that changes the basis vectors and therefore rotates the coordinate system cw through an angle φ.
Let φ = ωt, then U(R) rotates the coordinate system cw with angular frequency ω.  The coordinate system thus rotates with the magnetic field cw about the z-axis.

Consider an infinitesimal rotation and an infinitesimal time interval dt.
U(R)|ψ(dt)> = (I - (i/ħ)Szωdt)|ψ(dt))> = (I - (i/ħ)Szωdt)(I - (i/ħ)Hdt)|ψ(0)>.
U(R)|ψ(dt)> = (I - (i/ħ)Szωdt)(I - (i/ħ)(-ω0Sz - ω1(cosωt Sx - sinωt Sy))dt)|ψ(0)>.
U(R)|ψ(dt)> = (I - (i/ħ)((ω - ω0)Sz - ω1Sxcosωt + ω1Sysinωt)dt)|ψ(0)>.

S is a vector observable, its components transform under a rotation R as V’ = R-1V.
For a cw rotation R-1(φ) we have Vx’ = Vxcosφ - Vysinφ.
Therefore Sx’ = Sxcosωt - Sysinωt is the x-component of the observable S in a frame that rotates cw with angular frequency ω about the z-axis, i.e. in a frame that rotates with the magnetic field.
We therefore have  U(R)|ψ(dt)> = |ψ'(dt)> = (I - (i/ħ)((ω - ω0)Sz - ω1Sx’)dt)|ψ(0)>.
|ψ'(dt)> = (I  -(i/ħ)((ω - ω0)Sz’ - ω1Sx’)dt)|ψ’(0)>
|ψ'> denotes the state vector in a frame rotating cw about the z-axis with angular frequency ω.
In this frame the evolution operator is U(t,0) = exp(-(i/ħ)((ω - ω0)Sz’ - ω1Sx’)t).
Let ω = ω0
Then |ψ'(dt)> = (I + (i/ħ)ω1Sx’dt)|ψ’(0)>.  The infinitesimal evolution operator is I + (i/ħ)ω1Sx’dt and the evolution operator is U(t,0) = exp((i/ħ)ω1Sx’t).
The Hamiltonian therefore is H’ = -ω1Sx’.
The matrix of the Hamiltonian is

in the {|+>, |->} basis.

At t = 0 the lab frame and the frame that rotates with B coincide.  The density matrix at t = 0 in the {|+>, |->} basis is given by
.
Since H’ = -ω1Sx’ it is convenient to change the basis.  In the {|+>x, |->x} eigenbasis of Sx the density matrix is given by Uρ(0)U, where U is the unitary transformation that transforms {|+>, |->} into {|+>x, |->x}.
The matrix of U is

In the {|+>x, |->x} basis we have
.
Here Δp = p+ - p-.
The matrix of Sx in the {|+>x, |->x} basis is

Let us now calculate ρ(t) in the reference frame rotating with B
(d/dt)ρ'(t) = (1/(iħ))[H',ρ'(t)].
The matrix of H' in the {|+>x, |->x} basis is

We therefore have

11/dt = 0,  dρ22/dt = 0,  dρ12/dt = iω1ρ12,  dρ21/dt = -iω1ρ21.
This yields  ρ11(t) = 1/2,  ρ22(t) = 1/2,  ρ12(t) = -exp(iω1t) Δp/2,  ρ21(t) = exp(-iω1t) Δp/2.
We can now calculate <S’x>, <S’y>, and <S’z> in the rotating frame.
<S’x(t) > =Tr(ρ'(t)S’x).
.
In the {|+>x, |->x} basis the matrix of Sy is and the matrix of Sz is  .
In the rotating reference frame we therefore have

In the rotating frame the magnetization rotates clockwise about the x’-axis.  There are times when the magnetization has only a z-component, and times when it entirely lies in the plane perpendicular to the z-direction.

What do we observe in the lab frame?
S
is a vector observable.  The expectation values of the components of a  of a vector observable behave under rotation as though the observable were a classical vector.  We have
<Si'> = ∑iRij<Sj>,  <Si> = ∑iRij-1<Sj'>.
R(φ) is a ccw rotation about the z-axis.
The matrix of R-1 is

We therefore have (with ω = ω0)
<Sx> = <S’x>cosω0t + <Sy’>sinω0t = (ħ/2)Δp sinω1t sinω0t,
<Sy> = -<S’x>sinω0t + <Sy’>cosω0t = (ħ/2)Δp sinω1t cosω0t,
<Sz> = <S’z> = (ħ/2)Δp cosω1t.
The macroscopic magnetization M = nγ<S> of the water sample changes from nγ(ħ/2)Δp to -nγ(ħ/2)Δp and back.  Since B1 << B0 and the eigenstate of the unperturbed Hamiltonian are not degenerate, the eigenstates of the perturbed Hamiltonian are approximately equal to |+> and |->, and the eigenvalues are approximately -(ħω0/2) and +(ħω0/2).  When the magnetization changes from nγ(ħ/2)Δp to -nγ(ħ/2)Δp, many protons make transitions from the |+> to the |-> state.  They have to absorb energy from the field.  When the magnetization changes back to nγ(ħ/2)Δp, energy is dumped back into the field.

NMR
A continuous wave experiment can be performed by sweeping the frequency of the rotating field B1 through the resonance frequency ω0 and by monitoring the power output of the power supply that supplies the current that produces B1.  Oscillations in the power output indicate that the frequency is sweeping through the resonance frequency.  If γ is known, ω0 = γB0 yields B0.  In this way the strength of an unknown magnetic field can be measured.  If B0 is known, ω0 = γB0 yields γ.  Instead of sweeping the frequency ω of the rotating field, a continuous wave experiment can also be performed by sweeping B0 while holding ω constant.

Experiments can also be performed with pulsed magnetic fields.  Assume that at time t1 = π/(2ω1) we turn off B1.  Then for t > t1 H = -ω0Sz.  In the lab frame we have
<Sx> = (ħ/2)Δp sinω0t1, <Sy> = (ħ/2)Δp cosω0t1, <Sz> = 0  at t = t1.

At t = t1 the magnetization lies in the x-y plane.  How does it evolve?
To find the density matrix ρ(t1) in the lab frame in the {|+>, |->} basis we can use
<Sx> = Tr(ρ,Sx), or (ħ/2)Δp sinω0t= (ħ/2)(ρ12 + ρ21) = ħ Re(ρ12),  Re(ρ12) = ½Δpsinω0t1,
<Sy> = Tr(ρ,Sy), or (ħ/2)Δp cosω0t1 = i(ħ/2)(ρ12 - ρ21) = -ħ Im(ρ12),  Im(ρ12) = -½Δpcosω0t1,
<Sz> = Tr(ρ,Sz), or 0 = (ħ/2)(ρ11 - ρ22),  ρ11 = ρ22 = ½, since the trace is 1.
.

For t > t1 we have dρ11/dt = 0,  dρ22/dt = 0,  dρ12/dt = iω0ρ12,  dρ21/dt = -iω0ρ21,
which yields

For t > t1 we therefore have
<Sx(t)> = Tr(ρ(t),Sx) = (ħ/2)Δp sinωt,
<Sy(t)> = Tr(ρ(t)Sy) = (ħ/2)Δp cosωt,
<Sz(t)> = 0.
The magnetization lies in the x-y plane and rotates cw about the z-axis.  As the magnetization rotate about the z-axis, it will induce a current in a pick-up coil located with its axis along the x-axis.  Plotting this current versus time yields a sinusoidal wave form.  This current will decay as a function of time due to dephasing.  This is called free-induction decay.  The spin-spin relaxation time T2 is the time with which the transverse magnetization decays.  Mxy = Mxy0e-t/T2.  Loss of transverse magnetization is due to molecular interactions and field inhomogeneity.

At equilibrium the net magnetization is M = (0,0,Mz).  As we have seen, it is possible to make Mz = 0.  The time constant with which Mz returns back to its equilibrium value is called the spin-lattice relaxation constant T1.  Mz = M0(1 - e-t/T1).  Spin lattice relaxation is the result of the nuclei transferring energy to the surrounding molecules as thermal energy.

In pure water T1 ~ T2 ~ 2s - 3s.  In biological materials T2 < T1.

Reference:

The Basics of MRI
Note the chapters on spin physics and NMR spectroscopy