The goal is to come up with a time-dependent solution to the beam splitter Hamiltonian for a single photon, using a Taylor series approach, which I'll try to show the steps for below:
BTW - I don't have a lot of time now, but I figured I'd start you off on how I did it, and later I'll come back and make it more complete (Assuming I did it right, feel free to correct me!).
First we start with the time evolution operator, and evolve our initial state:
$|\psi(t) \rangle = e^{-i H t}|\psi_0\rangle$
Next we write this exponential operator as a Taylor series:
$|\psi(t) \rangle = \left(\hat{I} + (-\frac{it}{\hbar})\hat{H}+ \frac{1}{2!}(-\frac{it}{\hbar})^2\hat{H}^2 +... \right)|\psi_0\rangle$
Which of course has a representation in series notation as:
$|\psi(t) \rangle = \sum_k \frac{1}{k!}(-\frac{it}{\hbar})^k \hat{H}^k |\psi_0\rangle$
Now we plug in the Hamiltonian of interest, the beam splitter Hamiltonian, $ \hat{H} = a^\dagger b + b^\dagger a$:
$|\psi(t) \rangle = \sum_k \frac{1}{k!}(-\frac{it}{\hbar})^k (a^\dagger b + b^\dagger a)^k |\psi_0\rangle$
Now here we're going to specify our initial condition, which is that our photon is in one of the two modes to start out with. Let's say our state is $|0\rangle_a|1\rangle_b = |0, 1\rangle$. (So nothing in mode a, and a photon in mode b). Writing this out:
$|\psi(t) \rangle = \sum_k \frac{1}{k!}(-\frac{it}{\hbar})^k (a^\dagger b + b^\dagger a)^k |0, 1\rangle$
Now we will "pull out" a k to operate it on $|0, 1\rangle$:
$|\psi(t) \rangle = \sum_k \frac{1}{k!}(-\frac{it}{\hbar})^k (a^\dagger b + b^\dagger a)^{(k-1)}(a^\dagger b + b^\dagger a) |0, 1\rangle$
Now we can work out easily how this operator operates on the initial state using, $a^\dagger b |0, 1\rangle = |0, 1\rangle$ and $b^\dagger a |0, 1\rangle = 0$. This simplifies to:
$|\psi(t) \rangle = \sum_k \frac{1}{k!}(-\frac{it}{\hbar})^k (a^\dagger b + b^\dagger a)^{(k-1)} |1, 0\rangle$
Now we see that we have flipped from $|0, |1\rangle$ to $|1, |0\rangle$! If we pull out another k, then we end up back in our origional state. Since we are working with an infinite sum of these operators $\left(\hat{I} + (-\frac{it}{\hbar})\hat{H}+ \frac{1}{2!}(-\frac{it}{\hbar})^2\hat{H}^2 +... \right)$, we can see that the even terms of the sum will end up leaving the state in $|0, 1\rangle$, while the odd terms will change it to $|1, 0\rangle$. When we rewrite our terms into even and odd sums, we see that we get two infinite series. One of them you can transform into Sin() and the other turns into Cos(). Finally you get a form like $|\psi\rangle = -i Sin(it/\hbar)|1,0\rangle + Cos(it/\hbar)|0, 1\rangle$