Solving for the time-evolution operator in a periodically driven system

Question

I am looking at the Hamiltonian

$$H(t)=\begin{pmatrix} 0 & e^{i\Omega t}\\ e^{-i\Omega t}& 0\end{pmatrix}$$

I am trying to solve for the unitary operator $U(t,0)=\mathcal{T}\exp(-i\int_0^t dt'\, H(t'))$, where I have transformed this into $$ U(t,0)\approx\prod_{n=0}^{N}\exp\left( -iH(ndt)dt \right) $$ To calculate it with Mathematica, I then implemented a similar code to the one found here: Solving a time-dependent Schrödinger equation. I changed it a little, so that my code reads

U[H_,  ti_,  tf_,  n_] := 
  Module[{dt = N[(tf - ti)/n], Value = ({{1, 0},{0, 1}})}, 
    Do[Value = Dot[MatrixExp[-I*H[t]*dt], Value], {t, ti, tf, dt}]; 
    Value]

I then evaluated

U[λ/2 {{0, Exp[I Ω t]}, {Exp[-I Ω t], 0}}, 0, 20, 100]

However, this just gives me a long product of matrix exponentials:

MatrixExp[(0. - 0.2 I) {{0, 1/2 E^((0. + 20. I) Ω) λ}, {1/2 E^((0. - 20. I) Ω) λ, 0}}[20.]].
MatrixExp[(0. - 0.2 I) {{0, 1/2 E^((0. + 19.8 I) Ω) λ}, {1/2 E^((0. - 19.8 I) Ω) λ, 0}}[19.8]] ...

I really want to write it in the form

$$ H(t)=P(t)e^{-iH_Ft} $$

where the function $P(t)$ has a period of $T=2pi/\Omega$. This way, I can simply take $t=T$ and take a matrix log and get $H_F$. How can I get my product in this form and thus find the value of $H_F$? Am I doing something wrong with the code, or am I just missing a step?

Answer 1

Not knowing exactly what you want to do with the time evolution operator, I will simply give my version of how to obtain it in matrix form without using the manual discretization procedure you attempted, because I see no reason for it. I'm using NIntegrate so the complete code is much shorter and faster:

Omega = 2 Pi;
hamiltonian = {{0, Exp[I Omega*t]}, {Exp[-I Omega*t], 0}};
Manipulate[
 Module[{ψ},
  Row[{"U(", tMax, ") = ", MatrixForm@Chop@Transpose@Through[
        (ψ /. Table[
            First@NDSolve[{I D[ψ[t], t] == hamiltonian.ψ[t], 
                 ψ[0] == ψinit},  
                 ψ, {t, 0,tMax}], 
            {ψinit, IdentityMatrix[2]}
            ])[tMax]]}]],
 {{tMax, 1, "Time"}, 0., 20}]

This shows the matrix for $U(t,0)$ in the canonical basis, given by IdentityMatrix[2]. To get this matrix, NDSolve is applied twice: once for each unit vector. To display it, I vary the end time tMax in a Manipulate, and wrap the result in MatrixForm. Instead of the latter, you can now add further manipulations, such as getting the eigenvalues etc.

Edit

I think what you really want is the Floquet matrix $U_F$, given by the time evolution operator at $t=T$, whose powers $U_F^n$ give you the time evolution operator $U(n T)$. Here is a way to get that matrix for a particular choice of driving frequency Omega:

Omega = 2 Pi;
tMax = 2 Pi/Omega;
Clear[t, ψ];
hamiltonian = {{0, Exp[I Omega*t]}, {Exp[-I Omega*t], 0}};
UF = Transpose[(ψ[t] /. 
     Table[First@
       NDSolve[{I D[ψ[t], t] == 
          hamiltonian.ψ[t], ψ[0] == ψinit}, ψ[
         t], {t, 0, tMax}], {ψinit, IdentityMatrix[2]}]) /. 
   t -> tMax]

(*
==> {{0.987963 - 0.147405 I, 
  2.67424*10^-9 - 0.0469203 I}, {-2.67424*10^-9 - 0.0469203 I, 
  0.987963 + 0.147405 I}}
*)