Core Classes and Functions
MBParams Class
class MBParams:
"""
Parameters for Maxwell-Bloch simulations
Args:
L (float): Medium length in meters (default: 0.03)
Nz (int): Number of spatial grid points (default: 120)
tmax (float): Maximum time in seconds (default: 6e-6)
Nt (int): Number of time points (default: 4000)
lam_p (float): Probe wavelength in meters (default: 780e-9)
lam_c (float): Control wavelength in meters (default: 780e-9)
N (float): Atomic density in m^-3 (default: 5e16)
mu (float): Dipole moment in C·m (default: 3.0e-29)
gamma_eg (float): Excited state decay rate (default: 2π×6 MHz)
gamma_sg (float): Ground state coherence decay (default: 2π×0.2 MHz)
Delta_p (float): Probe detuning from resonance (default: 0)
delta_c (float): Two-photon detuning (default: 0)
sigma_v (float): Doppler width in m/s (default: 180)
Nv (int): Number of velocity classes (default: 11)
control_geom (str): 'co' or 'counter' propagating (default: 'co')
"""
Pulse Generation Functions
def gaussian_pulse(t, t0, tau_fwhm, area=0.5*2*pi*1e6):
"""
Generate Gaussian pulse envelope
Args:
t: Time array
t0: Pulse center time
tau_fwhm: Full width at half maximum
area: Pulse area (integral of Rabi frequency)
Returns:
Array of Rabi frequencies vs time
"""
def control_ramp_storage(t, Omega_c0=2*pi*10e6, t0_on=2.0e-6,
ramp=0.8e-6, hold=0.8e-6):
"""
Generate control field for storage/retrieval sequence
Args:
t: Time array
Omega_c0: Peak control Rabi frequency
t0_on: Time to start ramping down
ramp: Ramp duration
hold: Storage time
Returns:
Control field Rabi frequency vs time
"""
Propagation Functions
def mb_propagate_rk4(params, Op_in, Oc, include_doppler=True):
"""
Maxwell-Bloch propagation using 4th-order Runge-Kutta
Args:
params: MBParams instance
Op_in: Input probe Rabi frequency vs time
Oc: Control field Rabi frequency vs time
include_doppler: Whether to include Doppler averaging
Returns:
t: Time array
z: Spatial array
Op: Probe field Op[z,t] after propagation
"""
def mb_propagate_euler(params, Op_in, Oc, include_doppler=True):
"""
Maxwell-Bloch propagation using Euler method
(Similar arguments and returns as mb_propagate_rk4)
"""
Analysis Functions
def analysis_euler_vs_rk4(outdir, Nt=3000, Nz=120, Nv=9, sigma_v=180.0):
"""
Compare Euler and RK4 solvers
Outputs:
- Fig_Euler_vs_RK4.png: Comparison plot
- mb_compare_euler_rk4.csv: Numerical data
"""
def analysis_retrieval_heatmap(outdir, ramps_us, gammas_kHz,
Nt=2500, Nz=100, include_doppler=False):
"""
Generate retrieval efficiency heatmap
Args:
ramps_us: List of ramp durations in microseconds
gammas_kHz: List of dephasing rates in kHz
Outputs:
- Fig_Retrieval_Heatmap.png: Heatmap visualization
- retrieval_efficiency_heatmap.csv: Efficiency matrix
"""
def analysis_co_vs_counter_multi_sigma(outdir, sigmas,
Nt=2500, Nz=100, Nv=9):
"""
Compare co- and counter-propagating geometries
Args:
sigmas: List of Doppler widths in m/s
Outputs:
- Fig_Co_vs_Counter_multi_sigma.png: Comparison plot
- co_counter_sigma_*.csv: Data for each Doppler width
"""
Three-Level Λ System
Energy Level Structure
The system consists of three atomic levels:
- |g⟩: Ground state (stable)
- |e⟩: Excited state (decay rate γeg)
- |s⟩: Storage state (metastable)
Hamiltonian (Rotating Wave Approximation)
$$H = -\hbar\Delta_p |e\rangle\langle e| - \hbar(\Delta_p - \delta_c)|s\rangle\langle s| + \frac{\hbar\Omega_p}{2}(|e\rangle\langle g| + |g\rangle\langle e|) + \frac{\hbar\Omega_c}{2}(|e\rangle\langle s| + |s\rangle\langle e|)$$
Optical Bloch Equations
The density matrix elements evolve according to:
$$\dot{\rho}_{eg} = -\left(\frac{\gamma_{eg}}{2} + i\Delta_p'\right)\rho_{eg} + i\frac{\Omega_p}{2}(\rho_{gg} - \rho_{ee}) + i\frac{\Omega_c}{2}\rho_{sg}$$
$$\dot{\rho}_{sg} = -\left(\gamma_{sg} + i\delta'\right)\rho_{sg} + i\frac{\Omega_c}{2}\rho_{eg} - i\frac{\Omega_p}{2}\rho_{se}$$
$$\dot{\rho}_{ee} = -\gamma_{eg}\rho_{ee} + i\frac{\Omega_p}{2}(\rho_{ge} - \rho_{eg}) + i\frac{\Omega_c}{2}(\rho_{se} - \rho_{es})$$
Doppler Effect
For atoms moving with velocity v, the detunings are modified:
$$\Delta_p' = \Delta_p - k_p v \quad \text{(probe detuning)}$$
$$\delta' = (\Delta_p - \delta_c) - (k_p \mp k_c)v \quad \text{(two-photon detuning)}$$
where ± depends on co-propagating (+) or counter-propagating (-) geometry.
EIT Condition
Perfect transparency occurs when:
$$\Delta_p = 0, \quad \delta_c = 0, \quad \gamma_{sg} \to 0$$
The EIT linewidth is approximately:
$$\Delta\omega_{EIT} \approx \frac{|\Omega_c|^2}{\gamma_{eg}} + \gamma_{sg}$$
Computational Approach
Spatial Discretization
The medium is divided into Nz slices of thickness Δz = L/Nz. The probe field propagates according to:
$$\Omega_p(z + \Delta z, t) = \Omega_p(z, t) + i\eta\Delta z \langle\rho_{eg}(z, t)\rangle_v$$
Time Integration Methods
Euler Method (1st Order)
y(t + dt) = y(t) + dt * f(y(t), t)
Pros: Simple, fast per step
Cons: Accumulates error, can be unstable for stiff systems
4th-Order Runge-Kutta (RK4)
k1 = f(y(t), t)
k2 = f(y(t) + dt/2 * k1, t + dt/2)
k3 = f(y(t) + dt/2 * k2, t + dt/2)
k4 = f(y(t) + dt * k3, t + dt)
y(t + dt) = y(t) + dt/6 * (k1 + 2*k2 + 2*k3 + k4)
Pros: Much more accurate, better stability
Cons: 4x more function evaluations per step
Doppler Averaging
The velocity distribution is sampled at Nv points:
$$v_j = -3\sigma_v + \frac{6\sigma_v j}{N_v-1}, \quad j = 0, 1, ..., N_v-1$$
$$w_j = \exp\left(-\frac{v_j^2}{2\sigma_v^2}\right) / \sum_k w_k$$
Stability Considerations
- CFL Condition: Δt ≪ Δz/c for wave propagation
- Oscillation Period: Δt ≪ 2π/max(|Ω|) for Rabi oscillations
- Decay Time: Δt ≪ 1/γ for exponential decay
- Limiter: Coherences are clipped to |ρ| ≤ 5 to prevent numerical overflow
Performance Optimization
# Vectorized operations for all velocity classes
for iz in range(Nz):
# Solve atomic equations for all velocities simultaneously
rho_eg_all = solve_atomic_vectorized(Op[iz,:], Oc, all_velocities)
# Weighted average for propagation
rho_eg_avg = np.sum(weights * rho_eg_all, axis=0)
# Update probe field
Op[iz+1,:] = Op[iz,:] + 1j * eta * dz * rho_eg_avg
Getting Started
Installation
# Clone the repository
git clone https://github.com/alovladi007/Full-Maxwell-Bloch-Simulation-for-EIT-Based-Quantum-Optics.git
cd Full-Maxwell-Bloch-Simulation-for-EIT-Based-Quantum-Optics
# Create virtual environment
python3 -m venv venv
source venv/bin/activate # On Windows: venv\Scripts\activate
# Install dependencies
pip install -r requirements.txt
Basic Usage
# Run all analyses with default parameters
python src/high_fidelity_qo_suite.py --outdir output --euler_rk4 --heatmap --cocounter
# Run with custom parameters
python src/high_fidelity_qo_suite.py \
--outdir my_results \
--euler_rk4 --er_Nt 5000 --er_Nz 200 \
--heatmap --hm_doppler \
--cocounter --cc_sigmas "50,100,200,400"
Parameter Tuning Guide
Grid Resolution
- Nt (time points): Increase for longer simulations or finer time resolution
- Nz (spatial points): Increase for better spatial resolution or longer media
- Nv (velocity classes): Increase for more accurate Doppler averaging (9-15 typical)
Physical Parameters
# Typical ranges for Rb-87 vapor cell experiments
N = 1e15 to 1e17 # Atomic density (m^-3)
gamma_eg = 2π × 6 MHz # D2 line natural linewidth
gamma_sg = 2π × (0.1 to 10) kHz # Ground state dephasing
sigma_v = 50 to 500 m/s # Doppler width (temperature dependent)
Omega_c = 2π × (1 to 50) MHz # Control Rabi frequency
L = 1 to 10 cm # Cell length
Common Use Cases
1. EIT Spectroscopy Simulation
# Scan probe detuning to see EIT window
for delta_p in np.linspace(-50e6, 50e6, 100):
params = MBParams(Delta_p=2*pi*delta_p)
# Run simulation and extract transmission
2. Quantum Memory Optimization
# Find optimal control field parameters
ramp_times = np.linspace(0.1e-6, 2e-6, 20)
storage_times = np.linspace(0.5e-6, 5e-6, 20)
# Run parameter sweep and analyze efficiency
3. Temperature Effects Study
# Vary Doppler width to simulate different temperatures
temperatures = [300, 400, 500, 600] # Kelvin
for T in temperatures:
sigma_v = np.sqrt(k_B * T / m_Rb87)
# Run simulation with temperature-dependent Doppler
Troubleshooting
Common Issues
- Numerical instability: Reduce time step (increase Nt) or use RK4 instead of Euler
- Unphysical results: Check that γsg < γeg and |Ω| < γeg
- Slow performance: Reduce Nv for faster Doppler averaging, or disable with include_doppler=False
- Memory errors: Reduce grid size or process in chunks
