propagator in transition method is not unitary if delta is not zero

[savowe]

the propagator matrix U is not generally time-reversable, meaning U(tau)*U(-tau) \neq 1 It is right now only 1, if the two-photon detuning delta is zero. This means that it proably has something to do with it. Maybe the formula is wrong? Example:

import numpy as np
import matplotlib.pyplot as plt
import aisim as ais

pos_params = {
     'mean_x': 0.0,
     'std_x' : 3.0e-3, # cloud radius in m
     'mean_y': 0.0,
     'std_y' : 3.0e-3, # cloud radius in m
     'mean_z': 0.0,
     'std_z' : 0,        # ignore z dimension, its not relevant here
}
vel_params = {
     'mean_vx': 0.0,
     'std_vx' : ais.convert.vel_from_temp(3.0e-6), # cloud velocity spread in m/s at tempearture of 3 uK
     'mean_vy': 0.0,
     'std_vy' : ais.convert.vel_from_temp(3.0e-6), # cloud velocity spread in m/s at tempearture of 3 uK
     'mean_vz': 0.0,
     'std_vz' : ais.convert.vel_from_temp(0*.2e-6),       
}

print(ais.__version__)

atoms = ais.create_random_ensemble_from_gaussian_distribution(pos_params, vel_params, int(1e4))

center_Rabi_freq = 2*np.pi*12.5e3
r_beam = 29.5e-3 / 2 # 1/e^2 beam radius in m
intensity_profile = ais.IntensityProfile(r_beam, center_Rabi_freq)
wave_vectors = ais.Wavevectors( k1 = 8.052945514e6, k2 = -8.052802263e6)

Omega_eff = intensity_profile.get_rabi_freq(atoms.position)

phase = 0

delta = wave_vectors.doppler_shift(atoms)

# calculate Rabi frequency
Omega_R = np.sqrt(Omega_eff**2 + delta**2)
# beginning of pulse t0
t0 = atoms.time

tau = 20e-6
# calculate matrix elements
sin_theta = Omega_eff/Omega_R
cos_theta = -delta/Omega_R

U_ee = np.cos(Omega_R * tau / 2) - 1j * \
    cos_theta * np.sin(Omega_R * tau / 2)
U_ee *= np.exp(-1j * delta * tau/2)

U_eg = np.exp(-1j * (delta*t0 + phase)) * -1j * \
    sin_theta * np.sin(Omega_R * tau / 2)
U_eg *= np.exp(-1j * delta * tau/2)

U_ge = np.exp(+1j * (delta*t0 + phase)) * -1j * \
    sin_theta * np.sin(Omega_R * tau / 2)
U_ge *= np.exp(1j * delta * tau/2)

U_gg = np.cos(Omega_R * tau / 2) + 1j * \
    cos_theta * np.sin(Omega_R * tau / 2)
U_gg *= np.exp(1j * delta * tau/2)

propagator = np.array([[U_ee, U_eg], [U_ge, U_gg]], dtype='complex')
propagator = np.transpose(propagator, (2, 0, 1))
prop1 = propagator

tau = -tau
# calculate matrix elements
sin_theta = Omega_eff/Omega_R
cos_theta = -delta/Omega_R

U_ee = np.cos(Omega_R * tau / 2) - 1j * \
    cos_theta * np.sin(Omega_R * tau / 2)
U_ee *= np.exp(-1j * delta * tau/2)

U_eg = np.exp(-1j * (delta*t0 + phase)) * -1j * \
    sin_theta * np.sin(Omega_R * tau / 2)
U_eg *= np.exp(-1j * delta * tau/2)

U_ge = np.exp(+1j * (delta*t0 + phase)) * -1j * \
    sin_theta * np.sin(Omega_R * tau / 2)
U_ge *= np.exp(1j * delta * tau/2)

U_gg = np.cos(Omega_R * tau / 2) + 1j * \
    cos_theta * np.sin(Omega_R * tau / 2)
U_gg *= np.exp(1j * delta * tau/2)

propagator = np.array([[U_ee, U_eg], [U_ge, U_gg]], dtype='complex')
propagator = np.transpose(propagator, (2, 0, 1))
prop2 = propagator

print(np.matmul(prop1,prop2))

[savowe]

resolved