Rabi splitting with TLS

Here, we demonstrate the socket-free TLS workflow using the maxwelllink.SingleModeSimulation electromagnetic solver. By resonantly coupling one classical cavity mode to a two-level system (TLS), we aim to monitor the cavity coordinate and verify the expected Rabi splitting in the frequency domain.

1. Defining Molecule

We first create a Molecule instance using the non-socket mode, i.e., we directly initialize the TLS within the Molecule class:

import numpy as np
import maxwelllink as mxl

frequency_au = 0.242
mu12 = 187

molecule = mxl.Molecule(
    driver="tls",
    driver_kwargs={
        "omega": frequency_au,
        "mu12": mu12,
        "orientation": 2,
        "pe_initial": 0e-3,
    }
)

[Init Molecule] Operating in non-socket mode, using driver: tls

2. Defining the single mode cavity

Then, we create a SingleModeSimulation instance which defines the parameters for a single harmonic oscillator. The pre-defined molecule is also attached to this class for coupled light-matter simulations.

This single-mode cavity obeys the following equations of motion:

\[\ddot{q}_{\rm c} = -\omega_{\rm c}^{2}\, q_c - \varepsilon \sum_i \mu_i - \gamma_c\, p_{\rm c} + D(t),\]

where the effective electric field of this cavity mode is

\[E(t) = -\varepsilon q_{\rm c}(t) .\]

Here, \(\varepsilon = \frac{\omega_{\rm c}}{\sqrt{\epsilon_0 V}}\) is the coupling strength, and the sum runs over the selected molecular axis of all molecules attaced to the SingleModeSimulation class. All quantities are in atomic units. Dipole self-energy is excluded in the calculation of the effective electric field \(E(t)\).

coupling_strength = 5e-5
dt_au = 1e-1
damping_au = 0e-4
total_steps = 40960

sim = mxl.SingleModeSimulation(
    molecules=[molecule],
    frequency_au=frequency_au,
    coupling_strength=coupling_strength,
    damping_au=damping_au,
    coupling_axis="z",
    drive=0.0,
    dt_au=dt_au,
    qc_initial=[0, 0, 1e-5],
    record_history=True,
    # excluding dipole self-energy term for TLS model
    include_dse=False,
)

print(
    f"Configured SingleModeSimulation with omega_c = {frequency_au:.3f} a.u. "
    f"and g = {coupling_strength:.3f} a.u."
)

sim.run(steps=total_steps)

init TLSModel with dt = 0.100000 a.u., molecule ID = 0
Configured SingleModeSimulation with omega_c = 0.242 a.u. and g = 0.000 a.u.
[SingleModeCavity] Completed 1000/40960 [2.4%] steps, time/step: 1.21e-04 seconds, remaining time: 4.85 seconds.
[SingleModeCavity] Completed 2000/40960 [4.9%] steps, time/step: 1.13e-04 seconds, remaining time: 4.57 seconds.
[SingleModeCavity] Completed 3000/40960 [7.3%] steps, time/step: 1.16e-04 seconds, remaining time: 4.43 seconds.
[SingleModeCavity] Completed 4000/40960 [9.8%] steps, time/step: 1.06e-04 seconds, remaining time: 4.21 seconds.
[SingleModeCavity] Completed 5000/40960 [12.2%] steps, time/step: 9.50e-05 seconds, remaining time: 3.96 seconds.
[SingleModeCavity] Completed 6000/40960 [14.6%] steps, time/step: 1.11e-04 seconds, remaining time: 3.86 seconds.
[SingleModeCavity] Completed 7000/40960 [17.1%] steps, time/step: 1.08e-04 seconds, remaining time: 3.74 seconds.
[SingleModeCavity] Completed 8000/40960 [19.5%] steps, time/step: 1.07e-04 seconds, remaining time: 3.62 seconds.
[SingleModeCavity] Completed 9000/40960 [22.0%] steps, time/step: 8.67e-05 seconds, remaining time: 3.42 seconds.
[SingleModeCavity] Completed 10000/40960 [24.4%] steps, time/step: 1.09e-04 seconds, remaining time: 3.32 seconds.
[SingleModeCavity] Completed 11000/40960 [26.9%] steps, time/step: 1.10e-04 seconds, remaining time: 3.22 seconds.
[SingleModeCavity] Completed 12000/40960 [29.3%] steps, time/step: 1.10e-04 seconds, remaining time: 3.12 seconds.
[SingleModeCavity] Completed 13000/40960 [31.7%] steps, time/step: 1.17e-04 seconds, remaining time: 3.04 seconds.
[SingleModeCavity] Completed 14000/40960 [34.2%] steps, time/step: 1.24e-04 seconds, remaining time: 2.96 seconds.
[SingleModeCavity] Completed 15000/40960 [36.6%] steps, time/step: 1.22e-04 seconds, remaining time: 2.87 seconds.
[SingleModeCavity] Completed 16000/40960 [39.1%] steps, time/step: 1.34e-04 seconds, remaining time: 2.79 seconds.
[SingleModeCavity] Completed 17000/40960 [41.5%] steps, time/step: 1.51e-04 seconds, remaining time: 2.74 seconds.
[SingleModeCavity] Completed 18000/40960 [43.9%] steps, time/step: 1.14e-04 seconds, remaining time: 2.62 seconds.
[SingleModeCavity] Completed 19000/40960 [46.4%] steps, time/step: 1.11e-04 seconds, remaining time: 2.50 seconds.
[SingleModeCavity] Completed 20000/40960 [48.8%] steps, time/step: 1.13e-04 seconds, remaining time: 2.39 seconds.
[SingleModeCavity] Completed 21000/40960 [51.3%] steps, time/step: 1.11e-04 seconds, remaining time: 2.27 seconds.
[SingleModeCavity] Completed 22000/40960 [53.7%] steps, time/step: 1.13e-04 seconds, remaining time: 2.16 seconds.
[SingleModeCavity] Completed 23000/40960 [56.2%] steps, time/step: 1.05e-04 seconds, remaining time: 2.04 seconds.
[SingleModeCavity] Completed 24000/40960 [58.6%] steps, time/step: 1.13e-04 seconds, remaining time: 1.92 seconds.
[SingleModeCavity] Completed 25000/40960 [61.0%] steps, time/step: 1.05e-04 seconds, remaining time: 1.80 seconds.
[SingleModeCavity] Completed 26000/40960 [63.5%] steps, time/step: 1.13e-04 seconds, remaining time: 1.69 seconds.
[SingleModeCavity] Completed 27000/40960 [65.9%] steps, time/step: 8.78e-05 seconds, remaining time: 1.56 seconds.
[SingleModeCavity] Completed 28000/40960 [68.4%] steps, time/step: 1.05e-04 seconds, remaining time: 1.45 seconds.
[SingleModeCavity] Completed 29000/40960 [70.8%] steps, time/step: 1.09e-04 seconds, remaining time: 1.34 seconds.
[SingleModeCavity] Completed 30000/40960 [73.2%] steps, time/step: 9.06e-05 seconds, remaining time: 1.22 seconds.
[SingleModeCavity] Completed 31000/40960 [75.7%] steps, time/step: 8.95e-05 seconds, remaining time: 1.10 seconds.
[SingleModeCavity] Completed 32000/40960 [78.1%] steps, time/step: 9.68e-05 seconds, remaining time: 0.98 seconds.
[SingleModeCavity] Completed 33000/40960 [80.6%] steps, time/step: 9.84e-05 seconds, remaining time: 0.87 seconds.
[SingleModeCavity] Completed 34000/40960 [83.0%] steps, time/step: 1.06e-04 seconds, remaining time: 0.76 seconds.
[SingleModeCavity] Completed 35000/40960 [85.4%] steps, time/step: 9.12e-05 seconds, remaining time: 0.65 seconds.
[SingleModeCavity] Completed 36000/40960 [87.9%] steps, time/step: 1.08e-04 seconds, remaining time: 0.54 seconds.
[SingleModeCavity] Completed 37000/40960 [90.3%] steps, time/step: 1.09e-04 seconds, remaining time: 0.43 seconds.
[SingleModeCavity] Completed 38000/40960 [92.8%] steps, time/step: 1.13e-04 seconds, remaining time: 0.32 seconds.
[SingleModeCavity] Completed 39000/40960 [95.2%] steps, time/step: 1.16e-04 seconds, remaining time: 0.21 seconds.
[SingleModeCavity] Completed 40000/40960 [97.7%] steps, time/step: 1.16e-04 seconds, remaining time: 0.11 seconds.

3. Retrieve simulation observables

After the simulation, we can retrieve the TLS trajectory from molecule.additional_data_history together with the cavity coordinate q_c(t) stored directly by SingleModeSimulation.

population = np.array([entry["Pe"] for entry in molecule.additional_data_history])
tls_time_au = np.array([entry["time_au"] for entry in molecule.additional_data_history])

qc_history = np.array(sim.qc_history)[:, 2]  # z-component
energy_history = np.array(sim.energy_history)
time_history = np.array(sim.time_history)

print(
    f"Collected {population.size} TLS samples and {qc_history.size} cavity samples."
)

Collected 40960 TLS samples and 40960 cavity samples.

4. Inspect time-domain Rabi oscillations

Because the TLS is at resonance with the classical cavity mode, Rabi oscillations can be observed for this coupled system.

import matplotlib.pyplot as plt

plt.figure(figsize=(7, 4))
plt.plot(time_history, qc_history, label="Cavity coordinate $q_c(t)$")
plt.xlabel("time (a.u.)")
plt.ylabel("$q_c$ (a.u.)")
plt.title("Single-mode cavity response")
plt.legend()
plt.tight_layout()
plt.show()

plt.figure(figsize=(7, 4))
plt.plot(tls_time_au, population, label="TLS excited-state population")
plt.xlabel("time (a.u.)")
plt.ylabel("Pe")
plt.title("TLS population dynamics")
plt.legend()
plt.tight_layout()
plt.show()

plt.figure(figsize=(7, 4))
plt.plot(tls_time_au, energy_history, label="Total system (EM + TLS) energy")
plt.xlabel("time (a.u.)")
plt.ylabel("E (a.u.)")
plt.title("Total system energy dynamics")
plt.ylim(0, np.max(np.array(energy_history))*1.1)
plt.legend()
plt.tight_layout()
plt.show()

6. Fourier analysis of the cavity coordinate

Analytically, for this simple model system, we can calculate the assoicated Rabi splitting under the rotating wave approximation as:

\[\Delta_{\rm{Rabi}} = \varepsilon \mu_{12} \sqrt{\frac{2}{\hbar \omega }}\]

We can of course obtain the numerical Rabi splitting by Fourier transforming the cavity coordinate dynamics, which can be compared with the above analytical solution.

signal = qc_history - np.mean(qc_history)
if signal.size == 0:
    raise RuntimeError("No cavity data recorded; ensure the simulation was executed above.")

dt_sim = np.mean(np.diff(time_history)) if time_history.size > 1 else dt_au
fft_vals = np.fft.rfft(signal)
freqs = np.fft.rfftfreq(signal.size, d=dt_sim) * 2.0 * np.pi
spectrum = np.abs(fft_vals)

# analytical Rabi splitting
rabi_splitting = coupling_strength * mu12 * (2.0/frequency_au)**0.5

expected_peaks = np.array([
    frequency_au - rabi_splitting/2.0,
    frequency_au + rabi_splitting/2.0,
])

plt.figure(figsize=(7, 4))
plt.plot(freqs, spectrum, label="|FFT(q_c)|")
for idx, freq in enumerate(expected_peaks):
    label = r"$\omega_0 \pm \Delta/2$ (analytical)" if idx == 0 else None
    plt.axvline(freq, color="red", linestyle="--", label=label)
plt.xlabel("frequency (a.u.)")
plt.xlim(frequency_au*0.5, frequency_au*1.5)
plt.ylabel("spectral amplitude")
plt.title("Cavity-mode spectrum")
plt.legend()
plt.tight_layout()
plt.show()