PQLS: A High-Performance Python Library for Steady-State Simulation of Open Quantum Systems

I. INTRODUCTION

• Frequency-domain spectroscopy can require hundreds or thousands of independent steady-state solutions, and Doppler averaging evaluates responses over velocity distributions at every detuning.Sweeps over multiple physical parameters can make these repeated solves a dominant computational cost.Source: §I · p. 1
• Existing tools address complementary parts of quantum-system simulation, while constructing and solving independent configurations sequentially can still incur substantial computational overhead.The introduction names QuTiP, QuTiP-JAX, RydIQule, and ARC as established tools.Source: §I · pp. 1–2
• PQLS batches independent parameter configurations and uses JAX/XLA vectorization on CPUs and supported GPUs to reduce Python-level overhead and improve resource utilization.This approach replaces repeated sequential solver calls with a compiled batched computation.Source: §I · p. 2
• The three-tier API spans physical atomic ladder descriptions, intermediate Rabi-frequency, detuning, decay-rate, and topology inputs, and direct Hamiltonian and collapse-operator specification.Users can choose between physics-oriented modeling and direct numerical control.Source: §I · p. 2
• The high-level API uses ARC to determine atomic parameters, including transition frequencies, dipole matrix elements, Rabi frequencies, and decay rates.This is part of the physics-oriented atomic-ladder interface.Source: §I · p. 2
• Benchmarks against QuTiP, QuTiP-JAX, and RydIQule report substantial speedups for batched steady-state parameter sweeps.The introduction identifies numerical validation as well as performance benchmarking.Source: §I · p. 2

• Open systems are modeled with the Lindblad master equation, which describes density-matrix evolution under environmental interactions.The environment can be represented as a second system coupled to the system of interest.Source: §II · p. 2
• The Liouvillian superoperator combines the Hamiltonian commutator and dissipative contribution into one representation of system dynamics.The Hamiltonian encodes coherent interactions, while collapse operators represent relaxation processes such as spontaneous emission and blackbody-radiation-induced transitions.Source: §II · p. 2
• Under the rotating-wave approximation, the Hamiltonian for a coherently driven N-level system is an N × N Hermitian matrix with detunings on the diagonal and Rabi-frequency couplings off diagonal.Uncoupled states have zero coupling, Ωij = 0.Source: §II · p. 2

• PQLS uses JAX kernels optimized by XLA to minimize steady-state Lindblad-solver overhead and outperform NumPy or iterative QuTiP-based implementations.The paper points to Fig. 1 for an overview of the algorithm.Source: §III · p. 3

• The intermediate solve_custom_network API accepts explicitly defined transition and decay tuples, allowing non-ladder topologies such as networks with branching decays.It bypasses ARC and uses raw Rabi frequencies and detunings.Source: A. Physics-Aware Public API · pp. 3–4
• The high-level solve_quantum_ladder API supports N-level atomic ladder systems and obtains atomic properties and temperature-dependent decay rates through ARC.Users provide atomic topology and quantum numbers; supported atoms include 85Rb, 87Rb, and Cs.Source: A. Physics-Aware Public API · p. 3
• The low-level solve_steady_state API lets users directly provide batched Hamiltonian tensors and collapse operators.This level allows additional decay or dephasing mechanisms to be specified directly.Source: A. Physics-Aware Public API · pp. 3–4

IV. PQLS IMPLEMENTATION

• PQLS maps physical parameters to Hamiltonian and collapse-operator tensors, then uses JAX for its core execution.The solver reformulates the Lindblad dynamics through the Liouvillian to obtain the steady-state density matrix.Source: B. JAX Backend · p. 4, §IV · p. 4

B. JAX Backend

• The same PQLS code can run on CPUs or supported GPUs, and JAX enables differentiation of solver outputs with respect to model parameters.The paper identifies Fisher-information calculations and gradient-based driving-field optimization as possible applications of differentiability.Source: B. JAX Backend · p. 4

A. Superoperator Computation

• The steady-state condition becomes a linear system in Liouville space after vectorizing the N × N density matrix into a vector of size N^2.PQLS constructs the Liouvillian from the Hamiltonian and collapse operators using Kronecker products, replaces its first row with the trace constraint, and solves the resulting nonsingular system using JAX's LU-based linear solver.Source: A. Superoperator Computation · p. 4

B. Hardware-Accelerated Vectorization

• PQLS uses jax.vmap to process a batch of Hamiltonians represented as a tensor of shape (B, N, N), rather than invoking a solver separately for each parameter value.XLA chunks the independent systems across available CPU or GPU cores, enabling thousands of steady states in one compiled execution and reducing Python interpreter overhead.Source: B. Hardware-Accelerated Vectorization · pp. 4–5

A. Multi-Tiered API Validation

• The CI suite checks the high-level API on a 4-level 85Rb Rydberg RF ladder and the mid-level API on branching 3-level systems against sequential QuTiP solves.The low-level API is also tested on randomized Hermitian Hamiltonians for 2–6 levels and batches up to 128 configurations.Source: A. Multi-Tiered API Validation · p. 5
• 10^-11 is the maximum absolute difference allowed between the PQLS JAX backend and QuTiP across randomized trials.The suite separately sets a 10^-14 threshold for physical invariants, including trace conservation and Hermiticity.Source: A. Multi-Tiered API Validation · p. 5
• 10^-13 is the allowed negative eigenvalue floor for density-matrix positivity checks, accounting for floating-point round-off.The test suite requires eigenvalues to be at least −10^-13.Source: A. Multi-Tiered API Validation · p. 5
• 10^-6 is the minimum off-diagonal coherence magnitude required in validation cases involving coherent optical response.This threshold ensures the Rydberg EIT/ATS tests include meaningful coherences.Source: A. Multi-Tiered API Validation · p. 5

VI. SIMULATION EXAMPLES

• The examples use distinct API tiers: a physical ladder-system description for the 4-level 85Rb case and an intermediate API for the Λ-system with branching decay channels.The Λ-system connects two lower states to a common excited state.Source: A. AT Splitting · p. 5, B. Coherent Population Trapping (CPT) · p. 5

A. AT Splitting

• Both the high-level solve_quantum_ladder API and low-level solve_steady_state API produce the same numerical spectrum, including the ATS peaks.The low-level route accepts batched Hamiltonian tensors and spontaneous-decay jump operators.Source: A. AT Splitting · p. 5, Figure 3 · p. 6

B. Coherent Population Trapping (CPT)

• At zero two-photon detuning, destructive interference in the 3-level 87Rb Λ-system forms a dark state that suppresses probe absorption and creates a transparency window.Away from resonance, the dark-state condition is lost and absorption increases.Source: B. Coherent Population Trapping (CPT) · p. 5, Figure 4 · p. 6

VII. BENCHMARKS

• The 1,000-point benchmark records wall-clock execution time for a 4-level 87Rb steady-state simulation across parameter-space sizes.The reported workstation included an AMD Ryzen 9 7950X CPU and an NVIDIA GeForce RTX 4090 GPU.Source: A. Comparison with QuTiP and QuTiP-JAX · p. 6, §VII · p. 6

A. Comparison with QuTiP and QuTiP-JAX

• More than 2,700× speedup over both QuTiP and QuTiP-JAX was reported for a 1,000-point sweep, with agreement to reference solutions within the specified tolerance.PQLS compiles the sweep as one batched computation; Fig. 6 confirms numerical parity across varying RF field strengths.Source: A. Comparison with QuTiP and QuTiP-JAX · p. 6, Table I · p. 6

C. 2D Parameter Sweeps and Heatmaps

• Approximately 10 milliseconds on GPU sufficed for PQLS to compute 10,000 steady-state configurations and reproduce the 2D interference heatmaps for a 5-level 85Rb system.The sweep varies interference-field amplitude and coupling-laser detuning, reproducing the study from.Source: C. 2D Parameter Sweeps and Heatmaps · p. 6, Figure 5 · p. 7

B. Scaling Performance vs. RydIQule

• PQLS consistently outperforms RydIQule across the tested parameter counts, with its highest performance achieved using GPU acceleration.The comparison examines scalability from 10^3 to 10^6 parameters; the authors attribute the GPU result to saturating memory bandwidth.Source: B. Scaling Performance vs. RydIQule · p. 7
• At 1,000,000 systems, CPU time falls from 10.749 s with RydIQule to 3.271 s with PQLS; GPU PQLS takes 0.182 s, a reported speedup of up to 58×.The scaling benchmark also reports peak throughput above 5.4 million systems/s on GPU.Source: §VIII · p. 7, B. Scaling Performance vs. RydIQule · p. 7, Table II · p. 7

VIII. CONCLUSIONS

• For a 1,000-point sweep, PQLS reached 253,787 systems/s and an approximately 2,743× speedup over serial QuTiP.The conclusion also reports reduced CPU execution time versus RydIQule and a further reduction with GPU acceleration.Source: §VIII · p. 7

• The density matrix evolves by a Lindblad master equation with M decay channels, each represented by a collapse operator.Source: APPENDIX GENERAL CONSTRUCTION OF THE… · p. 8
• For a decay from state |sk⟩ to |dk⟩, the collapse operator is Ck = √γk|dk⟩⟨sk|, with γk the transition rate.Source: APPENDIX GENERAL CONSTRUCTION OF THE… · p. 8
• Row-wise vectorization expresses the master equation as a linear system using the Kronecker product and the N × N identity matrix.Source: APPENDIX GENERAL CONSTRUCTION OF THE… · p. 8
• The complete Liouvillian combines coherent and dissipative contributions.Source: APPENDIX GENERAL CONSTRUCTION OF THE… · p. 8