import numpy as np
from clerq import macro
# Rb87 D2 line
decay_rate = 6.065 # Gamma / 2π, MHz
wavelength = 780.24 # nmConverting from Experimental Units
clerq works in a consistent set of natural units, and the clerq.macro module converts the numbers you have in the lab (intensities, powers, temperatures, densities) into them.
- Frequencies are cyclic (\(\Omega/2\pi\), \(\Gamma/2\pi\), detunings), in MHz, and times are in µs.
- Lengths are in cm and number densities in cm\(^{-3}\).
- The
macrofunctions take wavelengths in nm, intensities in mW/cm\(^2\), powers in mW and beam waists in cm.
Every function returns a plain number ready to put in a system dictionary or a PropagationSolver.
Intensity to Rabi frequency
For a closed cycling transition the saturation intensity is \(I_\mathrm{sat} = \pi h c \Gamma / 3\lambda^3\) and the Rabi frequency is \(\Omega = \Gamma\sqrt{I / 2 I_\mathrm{sat}}\). Use factor to scale the dipole moment squared for a non-cycling transition, for example a squared Clebsch-Gordan coefficient.
i_sat = macro.saturation_intensity(decay_rate, wavelength) # mW/cm^2
peak_intensity = macro.gaussian_beam_peak_intensity(power=1.0, waist=0.1) # 1 mW, 1 mm waist
rabi_freq = macro.intensity_to_rabi_freq(peak_intensity, decay_rate, wavelength)
print(f"I_sat = {i_sat:.3f} mW/cm^2, peak = {peak_intensity:.1f} mW/cm^2, Omega/2π = {rabi_freq:.2f} MHz")I_sat = 1.669 mW/cm^2, peak = 63.7 mW/cm^2, Omega/2π = 26.49 MHz
macro.rabi_freq_to_intensity is the inverse. The result is the rabi_freq (or ampl of a time function) of a clerq Field.
Number density and interaction strength
The propagation equation is \(\partial_z\Omega = 2 i N g \rho_{01}\), where \(g/2\pi\) is macro.interaction_strength, per atom, and \(N\) is the number density. Multiply by the density to get the interaction_strengths entry. A weak resonant field is absorbed with optical depth \(N\sigma_0 L\), with \(\sigma_0 = 3\lambda^2/2\pi\).
n = 1e11 # cm^-3
length = 0.5 # cm
g = macro.interaction_strength(decay_rate, wavelength)
interaction_strengths = [n * g] # MHz / cm
od = macro.resonant_optical_depth(n, length, wavelength)
print(f"interaction_strengths = {interaction_strengths[0]:.3f} MHz/cm, OD = {od:.1f}")interaction_strengths = 440.726 MHz/cm, OD = 145.3
Doppler broadening
velocity_classes takes a thermal_width, the \(w\) in weights \(e^{-(\delta/w)^2}\). That is the most probable speed over the wavelength, which macro.thermal_width gives from the temperature and mass in amu. The FWHM is \(2\sqrt{\ln 2}\,w\).
w = macro.thermal_width(temperature=300.0, mass=86.909, wavelength=wavelength)
print(f"thermal_width = {w:.0f} MHz, FWHM = {2 * np.sqrt(np.log(2)) * w:.0f} MHz")thermal_width = 307 MHz, FWHM = 511 MHz
Pulse durations and areas
Time functions take times in µs and a cyclic amplitude in MHz. t_funcs.gaussian and t_funcs.sech accept either ampl or a pulse area n_pi (in units of \(\pi\)), so a 20 ns FWHM pulse is fwhm=0.02, and a pulse of known peak intensity uses the Rabi frequency from above as ampl.
Putting it together
A weak 1 µs pulse through a Rb87 D2 vapour (no Doppler broadening), with the field, optical depth and decay rate all taken from laboratory numbers. The Rabi frequency is far below \(\Gamma\), so the medium is in the linear regime and the resonant amplitude absorption of the solver should equal half the optical depth.
A Clebsch-Gordan coefficient scales the dipole moment squared once: either as factor in the macro functions, or in a Field’s factors, not both.
from clerq import PropagationSolver, spectral
n, length = 1.4e8, 1.0 # cm^-3, cm
weak_rabi_freq = macro.intensity_to_rabi_freq(0.01, decay_rate, wavelength) # 10 µW/cm^2
solver = PropagationSolver(
system={
"num_states": 2,
"decays": [{"channels": [[0, 1]], "rate": decay_rate}],
"fields": [
{
"coupled_levels": [[0, 1]],
"rabi_freq": weak_rabi_freq,
"rabi_freq_t_func": "gaussian",
"rabi_freq_t_args": {"ampl": 1.0, "fwhm": 1.0, "centre": 0.0},
}
],
},
t_min=-5.0,
t_max=15.0,
t_steps=800,
z_min=0.0,
z_max=length,
z_steps=10,
interaction_strengths=[n * macro.interaction_strength(decay_rate, wavelength)],
)
solver.solve()
freqs = spectral.freq_list(solver)
resonance = np.argmin(np.abs(freqs))
od = macro.resonant_optical_depth(n, length, wavelength)
print(f"solver: {spectral.absorption(solver, 0, -1)[resonance]:.4f}, OD / 2: {od / 2:.4f}")solver: 0.2029, OD / 2: 0.2035