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

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?

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],
ψ == ψinit},
ψ, {t, 0,tMax}],
{ψinit, IdentityMatrix}
])[tMax]]}]],
{{tMax, 1, "Time"}, 0., 20}] This shows the matrix for $U(t,0)$ in the canonical basis, given by IdentityMatrix. 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], ψ == ψinit}, ψ[
t], {t, 0, tMax}], {ψinit, IdentityMatrix}]) /.
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}}
*)