In many treatments in current journals you find authors using the interaction picture.  Assume the Hamiltonian of an arbitrary system is H₀(t), and the corresponding evolution operator is U₀(t,t₀).  We have
(iħ∂/∂t)U₀(t,t₀) = H₀(t)U₀(t,t₀),   U₀(t₀,t₀) = I.
Now assume that the system is perturbed in such a way that the Hamiltonian becomes
H(t) = H₀(t) + W(t).  For such a system the state vector in the interaction picture, |Ψ_I(t)> is defined from the state vector in the Schroedinger picture through
|Ψ_I(t)> = U₀^†(t,t₀)|Ψ_S(t)>.

How does |Ψ_I(t)> evolve?
iħ∂|Ψ_I(t)>/∂t = iħ∂/∂t(U₀^†(t,t₀)|Ψ_S(t)>)
=  iħ[∂/∂t(U₀^†(t,t₀)]|Ψ_S(t)>) + iħU₀^†(t,t₀)[∂/∂t(|Ψ_S(t)>]
=  -U₀^†(t,t₀)H_0s(t)|Ψ_S(t)>) + U₀^†(t,t₀)H_s(t)|Ψ_S(t)>.
Here we used
iħ∂/∂t U₀(t,t₀)  = H₀(t)U₀(t,t₀),  -iħ∂/∂t U₀^†(t,t₀) = U₀^†(t,t₀)H₀(t).

We now can write
iħ∂|Ψ_I(t)>/∂t = -U₀^†H_0S(t)U₀U₀^†|Ψ_S(t)> + U₀^†H_s(t)U₀U₀^†|Ψ_S(t)>
= -H_0I(t)Ψ_I(t)> + (H_0I(t) + W_I(t))|Ψ_I(t)> = W_I(t))|Ψ_I(t)>.
Therefore
iħ∂|Ψ_I(t)>/∂t = W_I(t))|Ψ_I(t)>.

We can rewrite this differential equation in the form of an integral equation.
d|Ψ_I(t)> = (iħ)^-1W_I(t))|Ψ_I(t)>dt.
∫_t0^t d|Ψ_I(t)>  = (iħ)^-1∫_t0^t W_I(t')|Ψ_I(t')>dt'.
|Ψ_I(t)> = |Ψ_I(t₀)> + (iħ)^-1∫_t0^t W_I(t')|Ψ_I(t')>dt'
   ¯¯_______________________________¯¯
This integral equation can be solved by iteration.
The ket |Ψ_I(t)> can be expanded in a power series of the form
|Ψ_I(t)> = {I + (iħ)^-1∫_t0^t dt'W_I(t') + (iħ)^-2∫_t0^t dt'W_I(t') ∫_t0^t' dt''W_I(t'') + ... }|Ψ₀(t)>.
The interaction picture assigns part of the time dependence to the state vectors, and part to the operators.  Following the method used to derive the equation for the evolution of the operators in the Heisenberg picture, we can derive the equation for the evolution of operators in the interaction picture.
dA_I/dt = (∂A_S/∂t)_I+  (iħ)^-1[A_I,H_0I].

When do we use the interaction picture?
The interaction picture is often used when H_0S is time-independent and W_S(t) is a small correction to H_0S. 
Assume that the problem governed by H_0S is already solved, either exactly, or in some approximation. 
Assume W_S(t) = 0 for t < 0.  Then |Ψ_I(0)> = |Ψ_S(0)>.  Neglecting higher order terms we have
|Ψ_I(t)> = |Ψ_I(t₀)> + (iħ)^-1∫_t0^t W_I(t')|Ψ_I(0>dt'
with W_I(t) = exp(iH₀t/ħ) W_S(t) exp(-iH₀t/ħ).

Let {|n>} be an orthonormal eigenbasis of H₀ and let the system be in the eigenstate |m> at t = 0, i.e. |Ψ_I(0)> = |m>.  We have H₀|m> = E_m|m>,
The probability P(E_k,t) of finding the system in the eigenstate |k> of H₀ at time t, i.e. the probability of measuring the eigenvalue E_k, is |<k|Ψ_I(t)>|^2..  (The predictions of quantum mechanics are independent of the representation.)
<k|Ψ_I(t)> = <k|m> + (iħ)^-1∫₀^t <k|W_I(t')|m>dt'
= <k|m> + (iħ)^-1∫₀^t exp(i(E_k - E_m)t/ħ) <k|W_S(t')|m>dt'.
With <k|m> = 0 we have
P(E_k,t) = (1/ħ²)|∫₀^t exp(i(E_k - E_m)t/ħ) <k|W_S(t')|m>dt'|².

This is the result of first order time dependent perturbation theory.  P(E_k,t) is the transition probability.

If W_S(t) is independent of time, i.e. if a constant small term is added to the Hamiltonian at t = 0 then
P(E_k,t) = (1/ħ²)|<k|W_S|m>|²| (exp(i(E_k - E_m)t/ħ - 1)/(i(E_k - E_m)/ħ)|²

= (1/ħ²)|<k|W_S|m>|²| (exp(-iE_mt/ħ - exp(-iE_kt/ħ)/(i(E_k - E_m)/ħ)|²
= (|<k|W_S|m>|²/(E_k - E_m)²) (2 sin(ω_kmt/2))²
= (|<k|W_S|m>|²/(E_k - E_m)²) 4 sin²(ω_kmt/2).
Here ω_km = (E_k - E_m)/ħ and we have used |e^iθ - e^iφ| = 2sin((θ - φ)/2).

Problem:

Consider a one-dimensional, infinitely deep well.  Let U(x) = 0, for 0 < x < a, and U(x) = ∞ everywhere else.  The eigenstates of H₀ = p²/(2m) + U(x) are Φ_n(x) = (2/a)^1/2sin(nπx/a) with eigenvalues E_n = n²π²ħ²/(2ma²).
Assume that at t = 0 you put a coin on the bottom of the well.
W(x) = W₀, a/4 < x < 3a/4 for t > 0, W(x) = 0 everywhere else, H = H₀ + W.
If at t = 0 the system is in the state Φ₃(x), what is the probability of finding it in Φ₁(x) at time t?

Solution:

• Concepts:
  The eigenfunctions of the infinite square well, time-dependent perturbation theory
• Reasoning:
  W is a small correction to H₀.  The problem governed by H₀ is already solved.
  W = 0  for  t < 0.  The probability P(E_k,t) of finding the system in the eigenstate |k> of H₀ at time t, i.e . the probability of measuring the eigenvalue E_k, is
  P(E_k,t) = (1/ħ²)|∫₀^texp((i/ħ)(E_k-E₃)t')<Φ_k|W|Φ₃>dt'|².
  P(E_k,t) is the transition probability.
  [Note: We are calculating the probability of finding the system in the ground state of the unperturbed Hamiltonian H₀, not of the perturbed Hamiltonian H.  We are calculating the probability that we find the system in the ground state after we take the coin out at time t.]
• Details of the calculation:
  P(E₁,t) = (1/ħ²)|∫₀^texp((i/ħ)(E₁-E₃)t')<Φ₁|W|Φ₃>dt'|²
  = |(<Φ₁|W|Φ₃>/ħ)∫₀^texp(iω₁₃t')dt'|²
  [∫₀^texp(iω₁₃t')dt' = (-i/ω₁₃)(exp(iω₁₃t) - 1)
  = (-i/ω₁₃)exp(iω₁₃t/2)(exp(iω₁₃t/2) - exp(-iω₁₃t/2))
  = (2/ω₁₃)exp(iω₁₃t/2)sin(ω₁₃t/2)]
  P(E₁,t) = (|<Φ₁|W|Φ₃>|²/(ħω₁₃)²) 4 sin²(ω₁₃t/2).
  ω₁₃ = 8π²ħ/(2ma²).
  <Φ₁|W|Φ₃> = (2W₀/a)∫_a/4^3a/4dx sin(πx/a)sin(3πx/a)
  = (2W₀/π)∫_π/4^3π/4 dx' sin(x')sin(3x')
  = (2W₀/π)[(sin(2x')/4 - sin(4x')/6))|_π/4^3π/4  = -W₀/π.
  
  P(E₁,t) = (2W₀/(ħπω₁₃))²sin²(ω₁₃t/2) = (m²a⁴W₀²/(4ħ⁴π⁶))sin²(2π²ħt/(ma²)).
  P(E₁,t) oscillates with time t.
  For t > 0 the Hamiltonian is H.  At t = 0 the system is not in an eigenstate of H but in a eigenstate of H₀, an operator that does not commute with H.  The probability of finding the system in an eigenstate of an operator that does not commute with the Hamiltonian changes with time.