How to solve this coupled linear ODE

Question

I've tried using mathematica and for some reason it's not working I think it's a pretty standard problem but for some reason I'm having a hard time This is essentially a variant of the Rabbi problem I'll really appreciate the help :) I will appreciate both explanation on how to solve this using mathematica or a manual solution \begin{align} i\hbar \frac{d\alpha}{dt}&=\frac{\hbar\omega_0}{2}\alpha+\frac{dE}{2}e^{i\omega t}\beta\\ i\hbar \frac{d\beta}{dt}&=-\frac{\hbar\omega_0}{2}\beta+\frac{dE}{2}e^{-i\omega t}\alpha \end{align}

Answer 1

I think you are trying to solve simultaneous linear differential equations of first order of the form $$\frac{dx}{dt}=ax+by \\ \frac{dy}{dt}=cx+dy.$$ I would highly recommend to go through the book "An Introduction to Dynamical system and chaos" (Springer) by G. C. Layek.

Answer 2

We can rewrite the system as $$ \frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}\alpha\\\beta\end{pmatrix} =\begin{pmatrix}-i\omega_0/2&-\frac{idE}{2\hbar}e^{i\omega t}\\-\frac{idE}{2\hbar}e^{-i\omega t}&i\omega_0/2\end{pmatrix} \begin{pmatrix}\alpha\\\beta\end{pmatrix} $$ and let us introduce $\tilde\alpha=e^{i\omega_0 t/2}\alpha$ and $\tilde\beta=e^{-i\omega_0 t/2}\beta$, giving $$ \frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}\tilde\alpha\\\tilde\beta\end{pmatrix} =\underbrace{-\frac{idE}{2\hbar}\begin{pmatrix}&-e^{i(\omega+\omega_0) t}\\e^{-i(\omega+\omega_0) t}&\end{pmatrix}}_{A(t)}\begin{pmatrix}\tilde\alpha\\\tilde\beta\end{pmatrix} $$ the $A$ here is much nicer to handle.

So we compute \begin{align*} [A]\text{-terms}&:\int_0^t A(\tau)\,\mathrm{d}\tau \\ &= \left(\frac{dE}{2\hbar(\omega+\omega_0)}\right)\begin{pmatrix}&e^{i(\omega+\omega_0)t}-1\\e^{-i(\omega+\omega_0)t}-1&\end{pmatrix}\\ [A^2]\text{-terms}&:\int_0^t A(\tau_1)\int_0^{\tau_1} A(\tau_2)\,\mathrm{d}\tau_2\,\mathrm{d}\tau_1\\ &= \left(\frac{dE}{2\hbar(\omega+\omega_0)}\right)^2 \begin{pmatrix}i(\omega+\omega_0)t+1-e^{i(\omega+\omega_0)t}\\&e^{-i(\omega+\omega_0) t}-1- i(\omega+\omega_0)t\end{pmatrix}\\ &\dots \end{align*} You should see that the $[A^n]$-term has two entries which are $\pm(dE/(2\hbar(\omega+\omega_0)))^n$ times the tail part of the series $e^{\pm i(\omega+\omega_0)t}$. Hence summing gives ultimately $$ I+\int_0^t A(\tau)\,\mathrm{d}\tau+\int_0^t A(\tau_1)\int_0^{\tau_1} A(\tau_2)\,\mathrm{d}\tau_2\,\mathrm{d}\tau_1+\dots = $$

$$ \begin{pmatrix} \frac1{2\omega_1} ( \omega_+ e^{i\omega_- t/2} - \omega_- e^{i\omega_+ t/2}) & \ast \\ \ast & \frac1{2\omega_1}( \omega_+ e^{i\omega_1 t/2} - \omega_- e^{-i\omega_1 t/2} ) \end{pmatrix} $$ where $\omega_1=\sqrt{(\omega+\omega_0)^2 + 4 (\frac{dE}{2\hbar})^2}$, $\omega_\pm= \omega+\omega_0 \pm \omega_1$ and I left two entries ($\ast$) for you to calculate.

Answer 3

I'll begin by recasting this problem in a more concretely physical way. Introducing $\Psi=(\alpha,\beta)^\top$, we have

$$i\hbar\frac{d\Psi}{dt}=i\hbar \frac{d}{dt}\begin{pmatrix} \alpha\\ \beta\end{pmatrix}=\begin{pmatrix}\hbar \omega_0/2 & e^{i\omega t}dE/2\\ e^{-i\omega t}dE/2 & -\hbar\omega_0/2\end{pmatrix}\begin{pmatrix} \alpha\\ \beta\end{pmatrix}=H(t)\Psi$$ which is the time-dependent Schrodinger equation with Hamiltonian $$\displaystyle H(t) = \begin{pmatrix}\hbar \omega_0/2 & e^{i\omega t}dE/2\\ e^{-i\omega t}dE/2 & -\hbar\omega_0/2\end{pmatrix}.$$ The Hamiltonian can be expressed in terms of Pauli matrices as

\begin{align} H(t)&=\frac{\hbar\omega_0}{2}\begin{pmatrix} 1& 0 \\ 0 & -1\end{pmatrix}+\frac{dE}{2}\cos\omega t\begin{pmatrix} 0 & 1 \\ 1 & 0\end{pmatrix}-\frac{dE}{2}\sin\omega t\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}\\ &=\frac{\hbar \omega_0}{2}\sigma_z+\frac{dE}{2}(\cos(\omega t)\sigma_x-\sin(\omega t)\sigma_y) \end{align} Note that the $dE$ portion amounts to rotating the operator $\sigma_x$ around the $z$-axis by an angle $-\omega$. Explicitly,

\begin{align} e^{i\omega t\sigma_z/2}\sigma_x e^{-i\omega t \sigma_z/2} &=\begin{pmatrix} e^{i\omega t/2} & 0 \\ 0 & e^{-i \omega t/2}\end{pmatrix} \begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix} \begin{pmatrix} e^{-i\omega t/2} & 0 \\ 0 & e^{i\omega t/2}\end{pmatrix}\\ &=\begin{pmatrix} 0 & e^{i\omega t} \\ e^{-i\omega t} & 0\end{pmatrix} \end{align} and therefore the Hamiltonian may be expressed as

$$H(t)=\frac{\hbar\omega_0}{2}+\frac{dE}{2}e^{i\omega t\sigma_z/2}\sigma_x e^{-i\omega t \sigma_z/2}$$ This amounts to an electron subject to an external magnetic field $\vec{B}=\frac12 (dE\cos\omega t,dE\sin\omega t,\hbar\omega_0)$ which rotates in some cone around the $z$-axis.

This suggests choosing a reference frame which rotates along with the magnetic field. To do this, we introduce $\Psi_r =e^{-i\omega t\sigma_z/2}\Psi$ ($r$ for rotated). Then the relevant Schrodinger equation is \begin{align} i\hbar \frac{d}{dt}\Psi_r &=i\hbar \frac{d}{dt}e^{-i\omega t\sigma_z/2}\Psi\\ &=\frac{\hbar\omega}{2} \sigma_z \Psi_r+e^{-i\omega t\sigma_z/2}i\hbar \frac{d}{dt}\Psi\\ &=\frac{\hbar\omega}{2} \sigma_z \Psi_r+e^{-i\omega t\sigma_z/2}H(t)\Psi\\ &=\frac{\hbar\omega}{2} \sigma_z \Psi_r+e^{-i\omega t\sigma_z/2}H(t)e^{i\omega t\sigma_z/2}\Psi_r\\ &=H_r\Psi_r \end{align} where we now have the time-independent Hamiltonian $\displaystyle H_r=\frac{1}{2}\hbar(\omega+\omega_0) \sigma_z+\frac{dE}{2}\sigma_x$ in the rotating frame. This Hamiltonian is solvable by hand. To facilitate this, let $$ \Omega=\frac12\sqrt{(dE/\hbar)^2+(\omega+\omega_0)^2},\quad \tan\phi = \frac{dE}{\hbar(\omega+\omega_0)},$$ with which we may write $H_r=\hbar\Omega [\cos(\phi)\sigma_z+\sin (\phi)\sigma_x].$ The same trick as before means that we may diagonalize the Hamiltonian as $$H_r = e^{-i\phi \sigma_y/2}\hbar \Omega \sigma_z e^{i\phi\sigma_y/2}.$$

We now solve the Schrodinger equation as $\Psi_r(t)=U_r(t)\Psi_r(t=0)$ where

\begin{align} U_r(t)=e^{-i H_r t/\hbar} &=e^{-i\phi \sigma_y/2}e^{-i \Omega \sigma_z t}e^{i\phi\sigma_y/2}\\ &=\begin{pmatrix} \cos\frac{\phi}{2} & -\sin \frac{\phi}{2}\\ \sin \frac{\phi}{2} & \cos\frac{\phi}{2}\end{pmatrix} \begin{pmatrix}e^{-i B_r/\hbar} & 0 \\ 0 & e^{i B_r/\hbar}\end{pmatrix} \begin{pmatrix} \cos\frac{\phi}{2} & \sin \frac{\phi}{2}\\ -\sin \frac{\phi}{2} & \cos\frac{\phi}{2}\end{pmatrix}\\ &=\begin{pmatrix} \cos \Omega t -i \cos\phi \sin\Omega t & -i \sin\phi \sin \Omega t\\ -i \sin\phi \sin\Omega t & \cos\Omega t+i \cos\phi \sin\Omega t \end{pmatrix} \end{align} Returning to the original reference frame, the solution takes the form

$$\Psi(t)=e^{i\omega t\sigma_z/2}\Psi_r(t)=e^{i\omega t\sigma_z/2}U_r(t)\Psi_r(0)=U(t)\Psi(0)$$ where the time-evolution operator is now \begin{align} U(t) &= e^{i\omega t\sigma_z/2}U_r(t)\\ &= \begin{pmatrix} e^{i\omega t/2} & 0 \\ 0 & e^{-i\omega t/2}\end{pmatrix} \begin{pmatrix} \cos \Omega t -i \cos\phi \sin\Omega t & -i \sin\phi \sin \Omega t\\ -i \sin\phi \sin\Omega t & \cos\Omega t+i \cos\phi \sin\Omega t \end{pmatrix}\\ &=\begin{pmatrix} e^{i\omega t/2}(\cos \Omega t -i \cos\phi \sin\Omega t) & -i e^{i\omega t/2}\sin\phi \sin \Omega t\\ -i e^{-i\omega t/2}\sin\phi \sin\Omega t & e^{-i\omega t/2}(\cos\Omega t+i \cos\phi \sin\Omega t) \end{pmatrix}\\ &=\begin{pmatrix} e^{i\omega t/2}(\cos \Omega t -i \Omega^{-1}(\omega+\omega_0) \sin\Omega t) & -i e^{i\omega t/2}\Omega^{-1}(dE/\hbar) \sin \Omega t\\ -i e^{-i\omega t/2}\Omega^{-1}(dE/\hbar) \sin\Omega t & e^{-i\omega t/2}(\cos\Omega t+i \Omega^{-1}(\omega+\omega_0) \sin\Omega t) \end{pmatrix} \end{align} This represents the full closed-form solution. As a sanity check, if $dE=0$ then $\Omega=\frac12( \omega+\omega_0)$ and thus

$$U(t)\to \begin{pmatrix} e^{i\omega t/2}e^{-i(\omega+\omega_0)t/2} & 0\\ 0 & e^{-i\omega t/2}e^{-i(\omega-\omega_0)t/2} \end{pmatrix} = \begin{pmatrix} e^{-i\omega_0 t/2} & 0\\ 0 & e^{-i\omega_0 t/2} \end{pmatrix} =e^{-i\omega_0 t \sigma_z/2} $$ which agrees with $H(t)= \dfrac{\hbar\omega_0}{2}\sigma_z$ in this limit.