#!/usr/bin/env python3 # SPDX-License-Identifier: MIT """Verified low-saturation Doppler-cooling and force benchmarks. The program computes seven related artifacts: 1. counterpropagating-beam force curves versus scaled velocity; 2. analytic and finite-difference friction coefficients versus detuning; 3. weak-field and representative shared-saturation Doppler temperatures; 4. force, friction, temperature, and capture proxies versus intensity; 5. a deterministic friction-diffusion relaxation benchmark; 6. an SI-scale rubidium-87 D2 example; 7. machine-readable validation and convention metadata. NumPy is the only dependency. The canonical detuning is atom minus laser, Delta = omega_0 - omega_L, so red detuning is positive. The optional finite-intensity curves use a declared shared-saturation denominator and are not a multilevel optical-molasses model. """ from __future__ import annotations import argparse import csv import json import math import platform from pathlib import Path from typing import Any, Callable import numpy as np HBAR = 1.054_571_817e-34 BOLTZMANN = 1.380_649e-23 RUBIDIUM87_MASS = 1.443_160_60e-25 RUBIDIUM87_D2_WAVELENGTH = 780.241_209_686e-9 RUBIDIUM87_D2_LINEWIDTH_HZ = 6.065e6 FORCE_DETUNINGS = (-0.5, 0.25, 0.5, 1.0, 2.0) SHARED_SATURATIONS = (0.01, 0.1, 0.5, 2.0) WEAK_SATURATION = 0.1 def value_slug(value: float) -> str: sign = "m" if value < 0.0 else "" text = f"{abs(value):g}".replace(".", "p") return sign + text def weak_force( velocity: np.ndarray | float, *, detuning: float, saturation: float, ) -> np.ndarray: """Return F/(hbar k Gamma) for two independent weak beams.""" velocity_array = np.asarray(velocity, dtype=float) plus = 1.0 + 4.0 * (detuning + velocity_array) ** 2 minus = 1.0 + 4.0 * (detuning - velocity_array) ** 2 return 0.5 * saturation * (1.0 / plus - 1.0 / minus) def weak_directional_forces( velocity: np.ndarray | float, *, detuning: float, saturation: float, ) -> tuple[np.ndarray, np.ndarray]: """Return signed +x and -x beam forces in the same normalization.""" velocity_array = np.asarray(velocity, dtype=float) plus = ( 0.5 * saturation / (1.0 + 4.0 * (detuning + velocity_array) ** 2) ) minus = ( -0.5 * saturation / (1.0 + 4.0 * (detuning - velocity_array) ** 2) ) return plus, minus def shared_saturation_force( velocity: np.ndarray | float, *, detuning: float, saturation_per_beam: float, ) -> np.ndarray: """Return a declared two-beam shared-saturation extension. The denominator uses s_total = 2 s_beam. This is a controlled toy extension for power broadening, not an exact multibeam Bloch solution. """ velocity_array = np.asarray(velocity, dtype=float) base = 1.0 + 2.0 * saturation_per_beam plus = base + 4.0 * (detuning + velocity_array) ** 2 minus = base + 4.0 * (detuning - velocity_array) ** 2 return 0.5 * saturation_per_beam * ( 1.0 / plus - 1.0 / minus ) def analytic_friction( detuning: np.ndarray | float, *, saturation_per_beam: float, shared_saturation: bool, ) -> np.ndarray: """Return alpha/(hbar k^2).""" detuning_array = np.asarray(detuning, dtype=float) base = 1.0 + ( 2.0 * saturation_per_beam if shared_saturation else 0.0 ) return ( 8.0 * saturation_per_beam * detuning_array / (base + 4.0 * detuning_array**2) ** 2 ) def numerical_friction( detuning: float, *, saturation_per_beam: float, shared_saturation: bool, step: float, ) -> float: """Return alpha/(hbar k^2) from a five-point derivative at v = 0.""" if step <= 0.0: raise ValueError("finite-difference step must be positive") force: Callable[[float], np.ndarray] if shared_saturation: force = lambda velocity: shared_saturation_force( velocity, detuning=detuning, saturation_per_beam=saturation_per_beam, ) else: force = lambda velocity: weak_force( velocity, detuning=detuning, saturation=saturation_per_beam, ) derivative = ( float(force(-2.0 * step)) - 8.0 * float(force(-step)) + 8.0 * float(force(step)) - float(force(2.0 * step)) ) / (12.0 * step) return -derivative def normalized_diffusion( detuning: np.ndarray | float, *, saturation_per_beam: float, shared_saturation: bool, ) -> np.ndarray: """Return D_p/(hbar^2 k^2 Gamma) for the declared 3D ledger.""" detuning_array = np.asarray(detuning, dtype=float) base = 1.0 + ( 2.0 * saturation_per_beam if shared_saturation else 0.0 ) return saturation_per_beam / ( base + 4.0 * detuning_array**2 ) def doppler_temperature_ratio( detuning: np.ndarray | float, *, saturation_per_beam: float, shared_saturation: bool, ) -> np.ndarray: """Return T/T_D, where k_B T_D = hbar Gamma / 2.""" detuning_array = np.asarray(detuning, dtype=float) if np.any(detuning_array <= 0.0): raise ValueError("temperature model requires positive red detuning") base = 1.0 + ( 2.0 * saturation_per_beam if shared_saturation else 0.0 ) return ( base + 4.0 * detuning_array**2 ) / (4.0 * detuning_array) def golden_section_minimum( function: Callable[[float], float], lower: float, upper: float, *, tolerance: float = 1.0e-13, ) -> tuple[float, float]: if not lower < upper: raise ValueError("optimizer interval must have lower < upper") inverse_phi = (math.sqrt(5.0) - 1.0) / 2.0 c = upper - inverse_phi * (upper - lower) d = lower + inverse_phi * (upper - lower) fc = function(c) fd = function(d) while upper - lower > tolerance: if fc <= fd: upper = d d = c fd = fc c = upper - inverse_phi * (upper - lower) fc = function(c) else: lower = c c = d fc = fd d = lower + inverse_phi * (upper - lower) fd = function(d) location = 0.5 * (lower + upper) return location, function(location) def rk4_relaxation( times: np.ndarray, *, initial_mean: float, initial_temperature_ratio: float, max_step: float, ) -> tuple[np.ndarray, np.ndarray]: """Integrate d mean/dtau=-mean and d theta/dtau=-2(theta-1).""" if max_step <= 0.0: raise ValueError("relaxation max step must be positive") states = np.empty((len(times), 2), dtype=float) states[0] = (initial_mean, initial_temperature_ratio) def rhs(state: np.ndarray) -> np.ndarray: return np.array( [-state[0], -2.0 * (state[1] - 1.0)], dtype=float, ) for index in range(1, len(times)): interval = float(times[index] - times[index - 1]) steps = max(1, math.ceil(interval / max_step)) step = interval / steps state = states[index - 1].copy() for _ in range(steps): k1 = rhs(state) k2 = rhs(state + 0.5 * step * k1) k3 = rhs(state + 0.5 * step * k2) k4 = rhs(state + step * k3) state = state + step * (k1 + 2*k2 + 2*k3 + k4) / 6.0 states[index] = state return states[:, 0], states[:, 1] def force_curve_rows( samples: int, ) -> tuple[list[dict[str, Any]], dict[str, Any]]: velocities = np.linspace(-3.0, 3.0, samples) rows: list[dict[str, Any]] = [] force_curves = { detuning: weak_force( velocities, detuning=detuning, saturation=WEAK_SATURATION, ) for detuning in FORCE_DETUNINGS } plus, minus = weak_directional_forces( velocities, detuning=0.5, saturation=WEAK_SATURATION, ) shared_curves = { saturation: shared_saturation_force( velocities, detuning=0.5, saturation_per_beam=saturation, ) for saturation in SHARED_SATURATIONS } for index, velocity in enumerate(velocities): row: dict[str, Any] = { "k_v_over_Gamma": velocity, "plus_beam_force_d0p5_s0p1": plus[index], "minus_beam_force_d0p5_s0p1": minus[index], "net_force_d0p5_s0p1": plus[index] + minus[index], } for detuning, values in force_curves.items(): row[ f"weak_force_d{value_slug(detuning)}_s0p1" ] = values[index] for saturation, values in shared_curves.items(): row[ "shared_force_d0p5_s" f"{value_slug(saturation)}" ] = values[index] rows.append(row) red = force_curves[0.5] blue = force_curves[-0.5] positive_mask = (velocities > 0.0) & (velocities <= 0.2) oddness = max( float(np.max(np.abs(values + values[::-1]))) for values in force_curves.values() ) cancellation = max( abs(float(values[samples // 2])) for values in force_curves.values() ) directional_sum_error = float( np.max(np.abs(plus + minus - red)) ) validation = { "max_force_oddness_error": oddness, "max_zero_velocity_force": cancellation, "max_directional_sum_error": directional_sum_error, "largest_red_small_velocity_force": float( np.max(red[positive_mask]) ), "smallest_blue_small_velocity_force": float( np.min(blue[positive_mask]) ), "checks": { "odd_force": oddness < 1.0e-15, "balanced_zero_force": cancellation < 1.0e-15, "directional_sum": directional_sum_error < 1.0e-15, "red_damps": bool(np.all(red[positive_mask] < 0.0)), "blue_antidamps": bool(np.all(blue[positive_mask] > 0.0)), }, } if not all(validation["checks"].values()): raise RuntimeError(f"force-curve validation failed: {validation}") return rows, validation def friction_rows( samples: int, ) -> tuple[list[dict[str, Any]], dict[str, Any]]: detunings = np.linspace(-2.0, 2.0, samples) weak_analytic = analytic_friction( detunings, saturation_per_beam=WEAK_SATURATION, shared_saturation=False, ) weak_numeric = np.array( [ numerical_friction( float(detuning), saturation_per_beam=WEAK_SATURATION, shared_saturation=False, step=0.001, ) for detuning in detunings ] ) shared = { saturation: analytic_friction( detunings, saturation_per_beam=saturation, shared_saturation=True, ) for saturation in SHARED_SATURATIONS } rows: list[dict[str, Any]] = [] for index, detuning in enumerate(detunings): row: dict[str, Any] = { "Delta_over_Gamma": detuning, "delta_laser_minus_atom_over_Gamma": -detuning, "alpha_weak_analytic_over_hbar_k2": weak_analytic[index], "alpha_weak_five_point_over_hbar_k2": weak_numeric[index], } for saturation, values in shared.items(): row[ "alpha_shared_s" f"{value_slug(saturation)}_over_hbar_k2" ] = values[index] rows.append(row) numeric_optimum, negative_maximum = golden_section_minimum( lambda detuning: -float( analytic_friction( detuning, saturation_per_beam=WEAK_SATURATION, shared_saturation=False, ) ), 1.0e-8, 2.0, ) analytic_optimum = 1.0 / (2.0 * math.sqrt(3.0)) derivative_error = float( np.max(np.abs(weak_numeric - weak_analytic)) ) validation = { "max_five_point_friction_error": derivative_error, "numeric_weak_friction_optimum_Delta_over_Gamma": ( numeric_optimum ), "analytic_weak_friction_optimum_Delta_over_Gamma": ( analytic_optimum ), "weak_friction_optimum_error": abs( numeric_optimum - analytic_optimum ), "maximum_weak_alpha_over_hbar_k2": -negative_maximum, "checks": { "finite_difference_friction": derivative_error < 2.0e-11, "friction_optimum": ( abs(numeric_optimum - analytic_optimum) < 2.0e-8 ), "red_positive_friction": bool( np.all(weak_analytic[detunings > 0.0] > 0.0) ), "blue_negative_friction": bool( np.all(weak_analytic[detunings < 0.0] < 0.0) ), }, } if not all(validation["checks"].values()): raise RuntimeError(f"friction validation failed: {validation}") return rows, validation def temperature_rows( samples: int, ) -> tuple[list[dict[str, Any]], dict[str, Any]]: detunings = np.linspace(0.02, 2.5, samples) weak_temperature = doppler_temperature_ratio( detunings, saturation_per_beam=WEAK_SATURATION, shared_saturation=False, ) weak_friction = analytic_friction( detunings, saturation_per_beam=WEAK_SATURATION, shared_saturation=False, ) weak_diffusion = normalized_diffusion( detunings, saturation_per_beam=WEAK_SATURATION, shared_saturation=False, ) shared_temperatures = { saturation: doppler_temperature_ratio( detunings, saturation_per_beam=saturation, shared_saturation=True, ) for saturation in SHARED_SATURATIONS } rows: list[dict[str, Any]] = [] for index, detuning in enumerate(detunings): row: dict[str, Any] = { "Delta_over_Gamma": detuning, "x_two_Delta_over_Gamma": 2.0 * detuning, "T_over_TD_weak": weak_temperature[index], "alpha_weak_over_hbar_k2": weak_friction[index], "Dp_weak_over_hbar2_k2_Gamma": weak_diffusion[index], "T_over_TD_from_Dp_over_alpha": ( 2.0 * weak_diffusion[index] / weak_friction[index] ), } for saturation, values in shared_temperatures.items(): row[ "T_over_TD_shared_s" f"{value_slug(saturation)}" ] = values[index] rows.append(row) numeric_optimum, minimum = golden_section_minimum( lambda detuning: float( doppler_temperature_ratio( detuning, saturation_per_beam=WEAK_SATURATION, shared_saturation=False, ) ), 0.01, 3.0, ) shared_optima: dict[str, Any] = {} largest_shared_location_error = 0.0 largest_shared_value_error = 0.0 for saturation in SHARED_SATURATIONS: location, value = golden_section_minimum( lambda detuning, s=saturation: float( doppler_temperature_ratio( detuning, saturation_per_beam=s, shared_saturation=True, ) ), 0.01, 4.0, ) expected_location = 0.5 * math.sqrt(1.0 + 2.0*saturation) expected_value = math.sqrt(1.0 + 2.0*saturation) largest_shared_location_error = max( largest_shared_location_error, abs(location - expected_location), ) largest_shared_value_error = max( largest_shared_value_error, abs(value - expected_value), ) shared_optima[f"s_{value_slug(saturation)}"] = { "numeric_Delta_over_Gamma": location, "analytic_Delta_over_Gamma": expected_location, "numeric_T_over_TD": value, "analytic_T_over_TD": expected_value, } ratio_identity_error = float( np.max( np.abs( weak_temperature - 2.0 * weak_diffusion / weak_friction ) ) ) validation = { "numeric_weak_temperature_optimum_Delta_over_Gamma": ( numeric_optimum ), "numeric_weak_minimum_T_over_TD": minimum, "weak_temperature_location_error": abs(numeric_optimum - 0.5), "weak_temperature_value_error": abs(minimum - 1.0), "largest_diffusion_friction_identity_error": ratio_identity_error, "largest_shared_optimum_location_error": ( largest_shared_location_error ), "largest_shared_optimum_value_error": largest_shared_value_error, "shared_optima": shared_optima, "checks": { "weak_Doppler_location": abs(numeric_optimum - 0.5) < 1.0e-8, "weak_Doppler_value": abs(minimum - 1.0) < 1.0e-13, "diffusion_friction_identity": ratio_identity_error < 1.0e-14, "shared_optimum_locations": ( largest_shared_location_error < 3.0e-8 ), "shared_optimum_values": ( largest_shared_value_error < 1.0e-13 ), }, } if not all(validation["checks"].values()): raise RuntimeError( f"temperature validation failed: {validation}" ) return rows, validation def intensity_rows( samples: int, ) -> tuple[list[dict[str, Any]], dict[str, Any]]: saturations = np.logspace(-3.0, 1.0, samples) rows: list[dict[str, Any]] = [] weak_fixed_temperature = float( doppler_temperature_ratio( 0.5, saturation_per_beam=WEAK_SATURATION, shared_saturation=False, ) ) for saturation in saturations: optimum_detuning = 0.5 * math.sqrt(1.0 + 2.0*saturation) minimum_temperature = math.sqrt(1.0 + 2.0*saturation) fixed_alpha = float( analytic_friction( 0.5, saturation_per_beam=float(saturation), shared_saturation=True, ) ) optimum_alpha = float( analytic_friction( optimum_detuning, saturation_per_beam=float(saturation), shared_saturation=True, ) ) fixed_temperature = float( doppler_temperature_ratio( 0.5, saturation_per_beam=float(saturation), shared_saturation=True, ) ) peak_velocity, negative_peak = golden_section_minimum( lambda velocity, s=float(saturation): float( shared_saturation_force( velocity, detuning=0.5, saturation_per_beam=s, ) ), 1.0e-10, 4.0, tolerance=2.0e-11, ) rows.append( { "saturation_per_beam": saturation, "total_saturation": 2.0 * saturation, "fixed_Delta_over_Gamma": 0.5, "alpha_fixed_over_hbar_k2": fixed_alpha, "T_fixed_over_TD": fixed_temperature, "temperature_optimal_Delta_over_Gamma": optimum_detuning, "minimum_T_over_TD": minimum_temperature, "alpha_at_temperature_optimum_over_hbar_k2": ( optimum_alpha ), "positive_peak_kv_over_Gamma_at_fixed_detuning": ( peak_velocity ), "peak_abs_force_over_hbar_k_Gamma": -negative_peak, } ) weak_intensity_independence_error = abs( weak_fixed_temperature - 1.0 ) shared_temperatures = np.array( [row["minimum_T_over_TD"] for row in rows] ) peak_forces = np.array( [row["peak_abs_force_over_hbar_k_Gamma"] for row in rows] ) validation = { "weak_temperature_at_Delta_half_over_TD": ( weak_fixed_temperature ), "weak_intensity_independence_error": ( weak_intensity_independence_error ), "lowest_shared_minimum_T_over_TD": float( shared_temperatures[0] ), "highest_shared_minimum_T_over_TD": float( shared_temperatures[-1] ), "lowest_peak_force": float(peak_forces[0]), "highest_peak_force": float(peak_forces[-1]), "checks": { "weak_temperature_normalization": ( weak_intensity_independence_error < 1.0e-15 ), "power_broadening_raises_minimum": bool( np.all(np.diff(shared_temperatures) > 0.0) ), "force_grows_over_scan": bool( peak_forces[-1] > peak_forces[0] ), "force_below_single_beam_ceiling": bool( np.all(peak_forces < 0.5) ), }, } if not all(validation["checks"].values()): raise RuntimeError(f"intensity validation failed: {validation}") return rows, validation def convergence_rows() -> tuple[list[dict[str, Any]], dict[str, Any]]: detuning = 0.5 saturation = 0.1 analytic = float( analytic_friction( detuning, saturation_per_beam=saturation, shared_saturation=False, ) ) steps = (0.08, 0.04, 0.02, 0.01, 0.005, 0.0025) values = [ numerical_friction( detuning, saturation_per_beam=saturation, shared_saturation=False, step=step, ) for step in steps ] errors = [abs(value - analytic) for value in values] ratios = [ errors[index] / errors[index + 1] for index in range(len(errors) - 1) ] rows = [ { "finite_difference_step_kv_over_Gamma": step, "alpha_numeric_over_hbar_k2": value, "alpha_analytic_over_hbar_k2": analytic, "absolute_error": error, "refinement_ratio_to_next": ( ratios[index] if index < len(ratios) else None ), } for index, (step, value, error) in enumerate( zip(steps, values, errors) ) ] validation = { "analytic_alpha_over_hbar_k2": analytic, "largest_error": max(errors), "smallest_error": min(errors), "minimum_refinement_ratio": min(ratios), "checks": { "error_decreases": all( upper < lower for lower, upper in zip(errors, errors[1:]) ), "fourth_order_regime": min(ratios) > 12.0, "finest_step_accuracy": errors[-1] < 1.0e-10, }, } if not all(validation["checks"].values()): raise RuntimeError( f"finite-difference convergence failed: {validation}" ) return rows, validation def relaxation_rows( samples: int, max_step: float, ) -> tuple[list[dict[str, Any]], dict[str, Any]]: times = np.linspace(0.0, 5.0, samples) mean_numeric, hot_numeric = rk4_relaxation( times, initial_mean=1.0, initial_temperature_ratio=10.0, max_step=max_step, ) _, cold_numeric = rk4_relaxation( times, initial_mean=0.0, initial_temperature_ratio=0.0, max_step=max_step, ) mean_exact = np.exp(-times) hot_exact = 1.0 + 9.0 * np.exp(-2.0 * times) cold_exact = 1.0 - np.exp(-2.0 * times) rows = [ { "time_over_velocity_damping_time": time, "mean_velocity_numeric_over_initial": mean_numeric[index], "mean_velocity_analytic_over_initial": mean_exact[index], "hot_temperature_numeric_over_equilibrium": ( hot_numeric[index] ), "hot_temperature_analytic_over_equilibrium": hot_exact[index], "cold_temperature_numeric_over_equilibrium": ( cold_numeric[index] ), "cold_temperature_analytic_over_equilibrium": ( cold_exact[index] ), } for index, time in enumerate(times) ] mean_error = float(np.max(np.abs(mean_numeric - mean_exact))) hot_error = float(np.max(np.abs(hot_numeric - hot_exact))) cold_error = float(np.max(np.abs(cold_numeric - cold_exact))) validation = { "max_mean_velocity_error": mean_error, "max_hot_temperature_error": hot_error, "max_cold_temperature_error": cold_error, "checks": { "mean_relaxation": mean_error < 2.0e-10, "hot_variance_relaxation": hot_error < 2.0e-9, "cold_variance_relaxation": cold_error < 2.0e-10, "hot_cools_monotonically": bool( np.all(np.diff(hot_numeric) < 0.0) ), "cold_heats_monotonically": bool( np.all(np.diff(cold_numeric) > 0.0) ), }, } if not all(validation["checks"].values()): raise RuntimeError(f"relaxation validation failed: {validation}") return rows, validation def rubidium_rows() -> tuple[list[dict[str, Any]], dict[str, Any]]: mass = RUBIDIUM87_MASS wavelength = RUBIDIUM87_D2_WAVELENGTH linewidth = 2.0 * math.pi * RUBIDIUM87_D2_LINEWIDTH_HZ wave_number = 2.0 * math.pi / wavelength saturation = 0.1 detuning = 0.5 alpha_dimensionless = float( analytic_friction( detuning, saturation_per_beam=saturation, shared_saturation=False, ) ) alpha = HBAR * wave_number**2 * alpha_dimensionless damping_time = mass / alpha doppler_temperature = HBAR * linewidth / (2.0 * BOLTZMANN) recoil_energy = HBAR**2 * wave_number**2 / (2.0 * mass) recoil_temperature = recoil_energy / BOLTZMANN recoil_velocity = HBAR * wave_number / mass resonant_velocity = detuning * linewidth / wave_number rms_velocity = math.sqrt(BOLTZMANN*doppler_temperature / mass) force_scale = HBAR * wave_number * linewidth peak_velocity, negative_peak = golden_section_minimum( lambda velocity: float( weak_force( velocity, detuning=detuning, saturation=saturation, ) ), 1.0e-10, 4.0, tolerance=2.0e-12, ) row = { "isotope": "rubidium-87", "transition": "D2 representative closed-two-level benchmark", "wavelength_m": wavelength, "linewidth_Hz": RUBIDIUM87_D2_LINEWIDTH_HZ, "linewidth_rad_per_s": linewidth, "mass_kg": mass, "wave_number_per_m": wave_number, "saturation_per_beam": saturation, "Delta_over_Gamma": detuning, "Doppler_temperature_K": doppler_temperature, "recoil_temperature_K": recoil_temperature, "Doppler_to_recoil_temperature_ratio": ( doppler_temperature / recoil_temperature ), "recoil_velocity_m_per_s": recoil_velocity, "Doppler_rms_velocity_m_per_s": rms_velocity, "resonant_velocity_m_per_s": resonant_velocity, "alpha_kg_per_s": alpha, "velocity_damping_time_s": damping_time, "single_beam_force_ceiling_N": 0.5 * force_scale, "pair_peak_abs_force_N": -negative_peak * force_scale, "pair_peak_positive_velocity_m_per_s": ( peak_velocity * linewidth / wave_number ), } validation = { "Doppler_temperature_microkelvin": ( 1.0e6 * doppler_temperature ), "recoil_temperature_microkelvin": ( 1.0e6 * recoil_temperature ), "velocity_damping_time_milliseconds": ( 1.0e3 * damping_time ), "resonant_velocity_m_per_s": resonant_velocity, "checks": { "Doppler_temperature_range": ( 145.0e-6 < doppler_temperature < 147.0e-6 ), "recoil_temperature_range": ( 0.17e-6 < recoil_temperature < 0.20e-6 ), "damping_time_range": ( 0.20e-3 < damping_time < 0.22e-3 ), "resonant_velocity_range": ( 2.3 < resonant_velocity < 2.5 ), }, } if not all(validation["checks"].values()): raise RuntimeError(f"rubidium validation failed: {validation}") return [row], validation def format_csv_value(value: Any) -> Any: if value is None: return "" if isinstance(value, (float, np.floating)): return f"{float(value):.16g}" if isinstance(value, (int, np.integer)): return int(value) return value def json_scalar(value: Any) -> Any: if isinstance(value, np.generic): return value.item() raise TypeError( f"Object of type {value.__class__.__name__} is not JSON serializable" ) def write_csv(path: Path, rows: list[dict[str, Any]]) -> None: if not rows: raise ValueError(f"cannot write empty CSV: {path}") with path.open("w", newline="", encoding="utf-8") as handle: writer = csv.DictWriter(handle, fieldnames=list(rows[0])) writer.writeheader() for row in rows: writer.writerow( { key: format_csv_value(value) for key, value in row.items() } ) def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--force-samples", type=int, default=1201) parser.add_argument("--friction-samples", type=int, default=801) parser.add_argument("--temperature-samples", type=int, default=1001) parser.add_argument("--intensity-samples", type=int, default=241) parser.add_argument("--relaxation-samples", type=int, default=501) parser.add_argument("--relaxation-max-step", type=float, default=0.005) parser.add_argument( "--output-dir", type=Path, default=Path("laser-cooling-simulation-output"), ) return parser.parse_args() def main() -> None: args = parse_args() for name in ( "force_samples", "friction_samples", "temperature_samples", "intensity_samples", "relaxation_samples", ): value = getattr(args, name) if value < 5: raise ValueError(f"{name.replace('_', '-')} must be at least 5") if value % 2 == 0: raise ValueError(f"{name.replace('_', '-')} must be odd") if args.relaxation_max_step <= 0.0: raise ValueError("relaxation-max-step must be positive") force, force_validation = force_curve_rows(args.force_samples) friction, friction_validation = friction_rows(args.friction_samples) temperature, temperature_validation = temperature_rows( args.temperature_samples ) intensity, intensity_validation = intensity_rows( args.intensity_samples ) convergence, convergence_validation = convergence_rows() relaxation, relaxation_validation = relaxation_rows( args.relaxation_samples, args.relaxation_max_step, ) rubidium, rubidium_validation = rubidium_rows() args.output_dir.mkdir(parents=True, exist_ok=True) paths = { "force": args.output_dir / "laser-cooling-force-curves.csv", "friction": args.output_dir / "laser-cooling-friction.csv", "temperature": args.output_dir / "laser-cooling-temperature.csv", "intensity": args.output_dir / "laser-cooling-intensity.csv", "convergence": args.output_dir / "laser-cooling-convergence.csv", "relaxation": args.output_dir / "laser-cooling-relaxation.csv", "rubidium": args.output_dir / "laser-cooling-rubidium.csv", "metadata": args.output_dir / "laser-cooling-metadata.json", } write_csv(paths["force"], force) write_csv(paths["friction"], friction) write_csv(paths["temperature"], temperature) write_csv(paths["intensity"], intensity) write_csv(paths["convergence"], convergence) write_csv(paths["relaxation"], relaxation) write_csv(paths["rubidium"], rubidium) metadata = { "scope": { "model": ( "closed two-level atom in weak optical molasses with a " "declared representative shared-saturation extension" ), "force_normalization": "F / (hbar k Gamma)", "velocity_coordinate": "u = k v / Gamma", "canonical_detuning": "Delta = omega_0 - omega_L", "red_detuning": "Delta > 0", "translated_detuning": "delta = omega_L - omega_0 = -Delta", "saturation": "s is the on-resonance value for one beam", "claim": ( "Doppler-force and friction-diffusion benchmark, not a " "multilevel or apparatus-specific cooling prediction" ), }, "weak_force": { "saturation_per_beam": WEAK_SATURATION, "detunings_Delta_over_Gamma": list(FORCE_DETUNINGS), "velocity_range_kv_over_Gamma": [-3.0, 3.0], "samples": args.force_samples, "validation": force_validation, }, "friction": { "detuning_range_Delta_over_Gamma": [-2.0, 2.0], "samples": args.friction_samples, "finite_difference": "five-point centered derivative", "validation": friction_validation, }, "temperature": { "definition": "k_B T = D_p / alpha", "Doppler_scale": "k_B T_D = hbar Gamma / 2", "diffusion_convention": "d Var(p_x) / dt = 2 D_p", "geometry": ( "three orthogonal weak beam pairs and isotropic " "spontaneous emission" ), "detuning_range_Delta_over_Gamma": [0.02, 2.5], "samples": args.temperature_samples, "validation": temperature_validation, }, "shared_saturation_extension": { "saturations_per_beam": list(SHARED_SATURATIONS), "total_saturation": "s_total = 2 s_beam", "denominator": "1 + s_total + 4 detuning^2", "status": "representative two-level extension", "intensity_samples": args.intensity_samples, "validation": intensity_validation, }, "finite_difference_convergence": convergence_validation, "relaxation": { "equations": [ "d mean_velocity / d tau = -mean_velocity", "d (T/T_eq) / d tau = -2 (T/T_eq - 1)", ], "time_range_tau": [0.0, 5.0], "samples": args.relaxation_samples, "integrator": "classical explicit Runge-Kutta order 4", "max_internal_step": args.relaxation_max_step, "validation": relaxation_validation, }, "rubidium_benchmark": { "isotope": "rubidium-87", "transition": "D2 representative closed-two-level benchmark", "validation": rubidium_validation, }, "validation": { "all_checks_passed": all( all(section["checks"].values()) for section in ( force_validation, friction_validation, temperature_validation, intensity_validation, convergence_validation, relaxation_validation, rubidium_validation, ) ), "sections": { "force": force_validation["checks"], "friction": friction_validation["checks"], "temperature": temperature_validation["checks"], "intensity": intensity_validation["checks"], "convergence": convergence_validation["checks"], "relaxation": relaxation_validation["checks"], "rubidium": rubidium_validation["checks"], }, }, "limitations": [ "closed two-level transition", "plane waves and balanced counterpropagating beams", "independent rates only in the weak-saturation benchmark", "local internal steady state at each velocity", "nonrelativistic classical center-of-mass motion", "dilute independent particles", "three-dimensional isotropic-emission diffusion ledger", "no hyperfine, Zeeman, polarization-gradient, or dark states", "no magnetic trap, spatial intensity profile, or gravity", "no branching, loss, reabsorption, collisions, or density effects", "no sub-Doppler, narrow-line, Raman, or sideband cooling", "no stochastic trajectory ensemble", ], "provenance": { "neutral_atom_cooling_proposal": { "authors": "T. W. Haensch and A. L. Schawlow", "title": "Cooling of gases by laser radiation", "journal": "Optics Communications 13, 68-69 (1975)", "doi": "10.1016/0030-4018(75)90159-5", }, "sub_Doppler_observation": { "authors": ( "P. D. Lett, R. N. Watts, C. I. Westbrook, " "W. D. Phillips, P. L. Gould, and H. J. Metcalf" ), "title": "Observation of atoms laser cooled below the Doppler limit", "journal": "Physical Review Letters 61, 169-172 (1988)", "doi": "10.1103/PhysRevLett.61.169", }, }, "runtime": { "python": platform.python_version(), "numpy": np.__version__, "platform": platform.platform(), "random_seed": None, }, "license": "MIT", "outputs": [path.name for path in paths.values()], } if not metadata["validation"]["all_checks_passed"]: raise RuntimeError("at least one validation check failed") paths["metadata"].write_text( json.dumps( metadata, default=json_scalar, indent=2, sort_keys=True, ) + "\n", encoding="utf-8", ) print( "force oddness error = " f"{force_validation['max_force_oddness_error']:.9e}" ) print( "friction derivative error = " f"{friction_validation['max_five_point_friction_error']:.9e}" ) print( "friction optimum Delta/Gamma= " f"{friction_validation['numeric_weak_friction_optimum_Delta_over_Gamma']:.9f}" ) print( "temperature optimum Delta/Gamma = " f"{temperature_validation['numeric_weak_temperature_optimum_Delta_over_Gamma']:.9f}" ) print( "rubidium Doppler temperature= " f"{rubidium_validation['Doppler_temperature_microkelvin']:.6f} microK" ) print( "relaxation RK4 max error = " f"{max(relaxation_validation['max_mean_velocity_error'], relaxation_validation['max_hot_temperature_error'], relaxation_validation['max_cold_temperature_error']):.9e}" ) print(f"outputs = {args.output_dir.resolve()}") print("validation = all checks passed") if __name__ == "__main__": main()