Speed up the evaluation of For loop

The time consuming step is not the For loop, but the many integrals. Therefore, calculate the indefinite integral first:

integ[t_] = Integrate[H0 // FullSimplify, t] 

Now you can calculate the table of integrals:

integs = Table[integ[time[[i + 1]]] - integ[time[[i]]], {i, 1, Length@time - 1}] 

With this we can get U0t;

U0t = MatrixExp[-I #] & /@ integs; 

And from this we can get the rest:

ρ1 = {{1/2, 0, I/2, 0}, {0, -(1/2), 0, I/2}, {-(I/2), 0, 1/2, 0}, {0, -(I/2), 0, -(1/2)}}; ρt = # . ρ1 . #\[ConjugateTranspose] & /@ U0t; Mxt = Tr[# . I1x] & /@ ρt; Myt = Tr[# . I1y] & /@ ρt; 

It takes already extremely long to compute integral. Note that your integrands are basically a functions times a constant $4 \times 4$ matrix. So just integrating the function and then multiplying by the matrix should help already. We can accelerate this further by employing the Fundamental Theorem of Calculus: We use Integrate to compute anti-derivatives symbolically and then compute differences of their function values to get the integrals. Finally we use TensorProduct to multiply the scala integrals with the matrix:

F[t_] = ComplexExpand[ Re[Integrate[\[Omega]D (\[Sqrt]2 Sin[2 \[Beta]] Cos[\[Omega]r t + \[Gamma]] - (Sin[\[Beta]])^2 Cos[2 \[Omega]r t + 2 \[Gamma]]), t]]]; G[t_] = ComplexExpand[ Re[Integrate[\[Omega]CS (-\[Sqrt]2 Sin[2 \[Beta]] Cos[\[Omega]r t + \[Gamma]] + (Sin[\[Beta]])^2 Cos[2 \[Omega]r t + 2 \[Gamma]]), t]]]; \[CapitalDelta]F = Differences[F[time]]; \[CapitalDelta]G = Differences[G[time]]; integrals = TensorProduct[\[CapitalDelta]F, (2 I1z . I2z - I1x . I2x - I1y . I2y)] + TensorProduct[\[CapitalDelta]G, (I1z + I2z)]; 

All entries of integrals compute in a fraction of a second.

And here some further performance improvements. Mathematica is an interpreted language, and those do typically not well with for loops. While you code would be nearly optimimal in C or C++, there are better ways to achive your goals in Mathematica. For example, you do not need to preallocate arrays for your results if you can use command Map (/@). Btw., Map is one of the reasons why I adopted Mathematica in the first place; it is so much easier to understand, less error prone and often faster at runtime than Do, For, or Table.

Next thing we can use is to let Dot operation thread over the list of U0t. Here we can use Compile with the option RuntimeAttributes -> {Listable}.

Finally, Tr[A.B] is always a bad idea, because A.B is an $O(n^3)$ operation while Dot[Flatten[A].Flatten[B]] does the same and only requires $O(n^2)$ time. By clever reshaping of \[Rho]t with Flatten, we can even compute all Dots with a single call to a matrix-vector Dot. The backend BLAS library is highly optimized and performs the task 35 times faster on my machine. (And I think the Flatten is actually for free here.)

\[Rho]1 = {{1/2, 0, I/2, 0}, {0, -(1/2), 0, I/2}, {-(I/2), 0, 1/2, 0}, {0, -(I/2), 0, -(1/2)}}; cf = With[{\[Rho]1 = \[Rho]1}, Compile[{{U, _Complex, 2}}, U . \[Rho]1 . U\[ConjugateTranspose], RuntimeAttributes -> {Listable}] ]; U0t = MatrixExp /@ (-I integrals); \[Rho]t = cf[U0t]; With[{\[Rho] = Flatten[\[Rho]t, {{1}, {2, 3}}]}, Mxt = \[Rho] . Flatten[I1x]; Myt = \[Rho] . Flatten[I1y]; ];