Normalization Constant in Time Evolution of Density Matrix Given the Hamiltonian:
$$%H = \omega \left(|0\rangle \langle1| + |1\rangle \langle0|  \right) = \begin{bmatrix}
0 & \omega \\ 
  \omega & 0 
\end{bmatrix}$$, I want to find the final state $\rho(t_f)$of the given density operator:
$$\rho(0) =|0\rangle \langle0| = \begin{bmatrix}
1 & 0 \\ 
 0 & 0 
\end{bmatrix} $$
To do so I started by stating:
\begin{equation}
\rho(t_f) = U\rho(0)U^\dagger\\
 U=e^{-i\frac{t_f}{\hbar}H} \approx 1-i\frac{t_f}{\hbar}H\; \; \Rightarrow \;\;U^\dagger \approx 1+i\frac{t_f}{\hbar}H 
\end{equation}
Although once I compute $\rho(t_f)$ using the above formula I obtain a non normalized state:
$$Tr(\rho(t_f))\neq 1 \; \; \forall \omega \neq 0$$
Of course this problem could be solved if out of nowhere I multiplied my $\rho(t_f)$ with a normalization constant N:
$$N = \frac{1}{Tr(\rho(t_f))}$$
My question is: is there something wrong with my thought process or calculations? Or do I really just have to introduce a new normalization constant? I would not mind an explanation in the option that the latter was the case(even if just as a reference).
I worked with it for a bit, and this is what I got:

P.S.
As suggested, I fully expand the U operator:
$$\%mathbf{U}=e^{-i\frac{t_f}{\hbar}\mathbf{H}} = \sum_n^\infty \left(\frac{c^n}{n!}\mathbf{H}^n \right)$$
Where for simplification I defined $c =i\frac{t_f}{\hbar}$.
By introducing a new operator denoted as $\mathbf{H}'$ ($\mathbf{H'} = \frac{1}{\omega}\mathbf{H}$), I notice the property:
$$\mathbf{H}^n=\left\{\begin{matrix}\omega^{n} \mathbf{I},& if \;\; n = even \\ 
\omega^n \mathbf{H}',& \; \; \; \; if \;\; n = odd
\end{matrix}\right.$$
Hence, the problem to solve becomes:
$$\mathbf{\rho}(t_f) = -\left(\sum_n^\infty \frac{c^n}{n!}\mathbf{H}^n \right)\rho(0) \left(\sum_n^\infty \frac{c^n}{n!}\mathbf{H}^n \right)$$
$$=%-\left( \sum_n^\infty\frac{c^{2n}}{2n!}\omega^{2n}\mathbf{I} + \frac{c^{2n+1}}{(2n+1)!}\omega^{2n+1}\mathbf{H'}\right)\rho(0)\left( \sum_n^\infty\frac{c^{2n}}{2n!}\omega^{2n}\mathbf{I} + \frac{c^{2n+1}}{(2n+1)!}\omega^{2n+1}\mathbf{H'}\right)$$
$$=-\left[ \sum_n^\infty \omega^{2n} \left( \frac{c^{2n}}{2n!}\mathbf{I} + \frac{c^{2n+1}}{(2n+1)!}\omega\mathbf{H'}\right )\right]\rho(0)\left[ \sum_n^\infty \omega^{2n} \left( \frac{c^{2n}}{2n!}\mathbf{I} + \frac{c^{2n+1}}{(2n+1)!}\omega\mathbf{H'}\right )\right]$$
$$=-\left( \sum_n^\infty \omega^{2n} \begin{bmatrix}
\frac{c^{2n}}{2n!} & \omega\frac{c^{2n+1}}{(2n+1)!}\\ 
 \omega\frac{c^{2n+1}}{(2n+1)!}&\frac{c^{2n}}{2n!} 
\end{bmatrix} \right)\left( \sum_n^\infty \omega^{2n} \begin{bmatrix}
\frac{c^{2n}}{2n!} & \omega\frac{c^{2n+1}}{(2n+1)!}\\ 
 0&0 
\end{bmatrix} \right)$$
$$%= \sum_n^\infty \omega^{4n}\begin{pmatrix}\frac{c^{4n}}{\left(2n\right)!\left(2n\right)!}&\omega\frac{c^{4n+1}}{\left(2n\right)!\left(2n+1\right)!}\\ \omega \frac{c^{4n+1}}{\left(2n\right)!\left(2n+1\right)!}& \omega^2\frac{c^{4n+2}}{\left(2n+1\right)!\left(2n+1\right)!}\end{pmatrix} = \mathbf{\rho}(t_f)$$
 A: If you expand $U$ to linear order in $t$, your density matrix will also only have trace one to linear order $t$, so $\mathrm{tr}(\rho(t))=1+O(t^2)$. As long as you get this, you did everything fine. Of course, your results will only be correct as long as the terms of order $t^2$ and higher will be small compared to the rest.
A: Additionally to Norbert Schuch's answer and comments, where he points out your problem of normalization, I want to add a brief note about the exact calculations that can be performed in this example.
First, for convenience, I want to take the factor $\omega$ out of the definition of the Hamiltonian. Now we have to notice (you also stated it in a comment) that
$$ H^{n}=\begin{cases} |0\rangle\langle 0| +  |1\rangle\langle 1| &\quad \text{for } n \text{ even} \\ H &\quad \text{for } n \text{ odd} \end{cases}  \quad .$$
To proceed, the exponential of an operator is defined by
$$e^{c\,H} \equiv \sum\limits_n \frac{c^n}{n!}\, H^n  \quad,$$
for $c\in\mathbb{C}$.
You can use this relation for your operators $U$ and $U^\dagger$. To make us of the elaborated properties of the Hamiltonian, you have to split the series into even and odd terms. You'll find a very simple expression for these operators. From this, you can calculate $U\, \rho(0) \,U^\dagger$ and thus find a form of $\rho(t)$.
Edit: Of course, if you calculate the exact $\rho(t)$, then $\mathrm{Tr}\rho(t) = 1$.
Edit 2: I think it would be easier to first consider only the expansion of $U$, which reads ($\hbar=1$):
$$U(t) = \underbrace{\sum\limits_{n=0}^{\infty} \frac{(-i\,\omega\,t)^{2n}}{(2n)!} H^{2n}}_{\text{even}} + \underbrace{\sum\limits_{n=0}^{\infty} \frac{(-i\,\omega\,t)^{2n+1}}{(2n+1)!} H^{2n+1}}_{\text{odd}}  \quad.$$
If you now use the properties of $H$ for the even and odd series, you will find
$$ U(t) = \sum\limits_{n=0}^{\infty} \frac{(-i\,\omega\,t)^{2n}}{(2n)!}\, \left(|0\rangle\langle 0| +|1\rangle\langle 1| \right)
 + \sum\limits_{n=0}^{\infty} \frac{(-i\,\omega\,t)^{2n+1}}{(2n+1)!}\, H   \quad.$$
Now you can simplify this expression with the help of the sine and cosine function, i.e. their respective series expansions. From there it is easy to obtain $U^{\dagger}(t)$ and also straightforward to calculate $\rho(t)$. Still, if you have questions, let me know.
