4. Oscylacje Rabiego

Na koniec rozważań nad modelem Jaynes’a-Cummings’a rozważymy jeszcze ogólną dynamikę modelu, rozwiązując równanie Schrödingera. Zobaczymy przy tym w jaki sposób łączone stany pola i atomu ewoluują oscylując pomiędzy stanami wzbudzen pola, a stanami wzbudzeń atomu. Oscylacje takie nazywane są oscylacjami Rabiego i są zjawiskiem bardzo ogólnym, znajdującym zastosowanie w bardzo wielu dziedzinach fizyki.

Jak już zauważyliśmy, Ewolucję można podzielić na podprzestrzenie rozpięte przez \(|{e,n}\rangle\), \(|{g,n+1}\rangle\), w których ewolucja w rzeczywistości jest po prostu dynamiką kubitu. Operator ewolucji w tej podprzestrzeni przyjmuje postać:
\begin{equation}
U_{n}(t) = \exp\Big(-\omega(n+\frac{1}{2})t\Big) \exp\Big( -i\hat{n}\cdot \vec{\sigma} (|\vec{n}|t) \Big)
\tag{21}
\end{equation}
gdzie \(\vec{n} = |\vec{n}|\hat{n}\) jest wektorem Blocha postaci:
\begin{equation}
\vec{n} = (\mathcal{G}\sqrt{n+1}, 0, \frac{1}{2}(\omega_{0}-\omega))^{T}
\tag{22}
\end{equation}
Ewolucja jest po prostu obrotem wektora Blocha na sferze Blocha wokół wektora \(\hat{n}\) z częstością kołową \(|\vec{n}|\). Możemy zapisać wyrażenie na częstość kołową obrotu wprost:
\begin{equation}
\Omega = \sqrt{\frac{1}{4}(\omega_{0}-\omega)^{2} + \mathcal{G}^{2}(n+1)}
\tag{23}
\end{equation}
Jeżeli jesteśmy w rezonansie, tzn. gdy \(\omega = \omega_{0}\), to częstość jest proporcjonalna do \(\sqrt{n+1}\). Oznacza to, że nawet jeżeli pole jest w stanie próżni (\(n= 0\)), to nadal ewolucja oscyluje pomiędzy stanem \(|{e,0}\rangle\) oraz \(|{g,1}\rangle\), czyli istnieje niezerowe prawdopodobieństwo emisji spontanicznej.

Nie będziemy w tym miejscu rozwiązywać równania Schrödingera analitycznie, zachęcamy czytelnika, żeby spróbował przeprowadzić obliczenia samodzielnie w ramach ćwiczenia. Skupimy się na numerycznej symulacji, dzięki której zobaczymy główne idee związane z oscylacjami Rabiego. Przed rozpoczęciem obliczeń spróbujmy przewidzieć wyniki. Chcielibyśmy otrzymać prawdopodobieństwa przejść między stanami \(|{e,n}\rangle\) oraz \(|{g,n+1}\rangle\). Ponieważ wektor Blocha dla operatora ewolucji, ma tylko składowe w kierunkach \(x, z\), natomiast rozpatrywana baza związana jest z kierunkiem \(z\), to przewidujemy, że cała ewolucja będzie oscylacjami prawdopodobieństw. Dodatkowo, amplituda tych oscylacji powinna być maksymalna w rezonansie, tzn. gdy \(\omega_{0} = \omega\), ponieważ wtedy nasz stan będzie obracał się dookoła osi \(x\).

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm
from matplotlib import cm

sx = np.array([
    [0,1],
    [1,0]
])

sy = np.array([
    [0,-1j],
    [1j,0]
])

sz = np.array([
    [1,0],
    [0,-1]
])

S = [sx,sy,sz]

def generate_U(n : int, delta : float, G : float, om : float, dt : float) -> np.ndarray:
    """
    Argumenty:
      - n : liczba fotonów okreslająca podprzestrzeń
      - delta : om_0 - om szerokość przerwy energetycznej
      - G : stała sprzężenia oddziaływania atomu z polem
      - om : częstośc pola
      - dt : krok czasowy wyznaczający dokładność symulacji
    """

    n_vec = [G * np.sqrt(n+1), 0 , delta / 2]

    # faza dynamiczna
    U1 = np.exp(-1j * om * (n+1/2) * dt)
    # właściwy operator ewolucji
    U2 = expm(- 1j * sum(S_ * n_ for S_, n_ in zip(n_vec, S)) * dt)

    return U1 * U2

def get_probabilities(U : np.ndarray, N : int) -> tuple:
    """
    Argumenty:
      - U : operator ewolucji
      - N : liczba kroków symulacji
    """

    psi_0 = np.array([1,0]).reshape(-1,1)

    # listy zawierająceprawdopodobieństwa dla stanów |e>|n> |g>|n+1>
    P_e = []
    P_g = []

    for _ in range(N):
      psi_0 = U @ psi_0
      P_e.append(np.abs(psi_0[0,0])**2)
      P_g.append(np.abs(psi_0[1,0])**2)

    return np.array(P_e), np.array(P_g)

n = 1
delta_s = [0. , 1 , 5]
G = 1.
om = 1
dt = 1e-2

t = np.arange(0, 5, dt)

fig = plt.figure(figsize = (8,3))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)

for i, delta in enumerate(delta_s):
U_ = generate_U(n, delta, G, om, dt)
P_e, P_g = get_probabilities(U_, len(t))
ax1.plot(t, P_e, color = cm.rainbow(i / len(delta_s)), lw = 2, alpha = 1, label = r"$\Delta$" + f"={delta}")
ax2.plot(t, P_g, color = cm.rainbow(i / len(delta_s)), lw = 2, alpha = 1, label = r"$\Delta$" + f"={delta}")

ax1.set_title(r"$|\langle e,n | \psi(t) \rangle|^{2}$")
ax1.set_xlabel("t")
ax1.set_ylabel("P")
ax2.set_title(r"$|\langle g,n+1 | \psi(t) \rangle|^{2}$")
ax2.set_xlabel("t")
ax2.legend(loc = "upper right")

plt.show()

Rysunek 4. Oscylacje Rabiego w modelu JC.

Przyglądając się numerycznym wynikom, możemy zaobserwować, iż rzeczywiście, oscylacje osiągają maksymalną amplitudę dla przypadku rezonansowego, natomiast zwiększanie odchylenia częstości pola \(\omega\) od częstości wyznaczonej przez szerokość przerwy energetycznej między stanami, amplituda drgań zimniejsza się i przejścia między stanami stają się mniej prawdopodobne.