QArray+: A physics-informed GPU-accelerated simulator for quantum dot arrays

QArray+: A physics-informed GPU-accelerated simulator for quantum dot arrays Pranav Vaidhyanathan,1, 2, ∗ Barnaby van Straaten,3, ∗, † Alice Petrillo,3 Rahul Marchand,2 Edwin De Nicolo,2 Menno Veldhorst,3 Brucek Khailany,1 Taylor L. Patti,1 and Natalia Ares2, ‡

arXiv:2609.02736v1 [cond-mat.mes-hall] 2 Sep 2026

2

1 NVIDIA, Santa Clara, California 95051, USA Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, United Kingdom 3 QuTech and Kavli Institute of Nanoscience, Delft University of Technology, P.O. Box 5046, 2600 GA Delft, The Netherlands (Dated: September 3, 2026)

Semiconductor quantum-dot arrays are a compelling platform for scalable quantum technologies, yet their practical operation is hindered by the complexity of tuning large-scale devices. Existing automation tools rely on simplified physical models—such as constant-capacitance approximations and equilibrium Hubbard models—which assume instantaneous relaxation to a steady state. These frameworks fail in experimentally critical regimes where measurement rates exceed tunneling dynamics, necessitating more sophisticated non-equilibrium control strategies. To bridge this gap, we introduce QArray+, an extension of the QArray framework that incorporates gate-dependent tunnel coupling and a quantum open-system description of dissipative processes. This approach enables the unified simulation of coherent interdot charge-state hybridization and the non-equilibrium latching dynamics essential for training robust machine-learning models for automated device operation. Implemented in JAX with GPU acceleration, QArray+ scales across GPUs and multi-node systems. For example, a charge stability diagram for a 100 × 100 grid of gate-voltages over 64 dots can be computed in ∼ 0.17 s on multiple GPUs. Since interdot interactions are short-ranged and the corresponding tuning corrections are local, simulations at these scales capture the physics relevant to even larger devices. These capabilities support high-throughput dataset generation for automated device tuning.

I.

INTRODUCTION

Semiconductor quantum dots (QDs) are a compelling platform for scalable quantum technologies, offering highfidelity qubit operations and compatibility with modern semiconductor fabrication [1–4]. The realization of highfidelity quantum circuits based on large QD arrays requires their tuning, control, and operation to be automated. To this end, machine-learning-based strategies for charge tuning have shown strong promise, but their development depends on realistic simulation environments and large, labeled datasets [5–14]. Since acquiring such datasets experimentally is impractical, fast physicsbased simulators have become an appealing method for creating automated tuning pipelines [8, 11, 15–20]. Existing tools, typically based on constant-capacitance models or simplified Hubbard models, have proven effective for generating charge-stability diagrams (CSDs) used in the training of ML models [6–8]. These simulators generate charge-stability diagrams by determining the system’s steady-state charge configuration, usually assumed to be the one that minimizes the state’s electrostatic energy at each point in gate-voltage space. However, this implicitly assumes that the dots’ charge configuration can relax to that steady state on timescales much faster than the measurement timescale. This assumption breaks down whenever charge relaxation is

∗ These authors contributed equally. † b.l.vanstraaten@tudelft.nl ‡ natalia.ares@eng.ox.ac.uk

slow on the measurement timescale. In particular, the coupling of a dot to its reservoir is exponentially suppressed with distance, so that relaxation times can become extremely long, leading to slow charge dynamics and hysteresis [1, 21–24]. Overcoming these limitations requires simulators capable of capturing a broader range of physical phenomena while remaining computationally efficient. Specifically, advanced frameworks must be able to simulate non-equilibrium latching, coherent interdot tunnel coupling for charge hybridization, and the simultaneous presence of both via an open-quantum-system approach. Furthermore, realistic modeling demands accommodating for effects such as non-linear electrostatics through voltage-dependent capacitances, supporting isolated modes where dots are entirely decoupled from reservoirs, and utilizing GPU acceleration to keep these complex, large-scale calculations tractable. We introduce QArray+, a successor to QArray [15]. By incorporating both interdot tunnel coupling and finite reservoir tunneling rates, QArray+ extends the original capacitance-based framework to capture a broader range of experimentally observable phenomena, such as coherent interdot hybridization and non-equilibrium charge dynamics. This physics-informed simulator is explicitly designed for acceleration and scalability. Because it is implemented in JAX with just-in-time compilation and GPU execution, it is easily executable in a single GPU as well as multi-GPU and multi-node environments. Furthermore, its electrostatic core uses a polynomial reformulation that reduces the per-point compute cost of generating charge-stability diagrams, ultimately enabling dataset-scale synthesis. The package provides three com-

2 Framework Feature

QDSim QArray QDFlow QDarts SimCATS RF-Squad QArray+

Simulate tunnel coupling Simulate latching Simulate latching and tunnel coupling Simulate voltage-dependent capacitances Simulate isolated mode Simulate RF-reflectometry readout Utilize GPU acceleration Utilize multi-GPU acceleration

× × × × × × × ×

× ◦ × × ✓ × ✓ ×

× ◦ × ✓ × × × ×

✓ × × ✓ × × × ×

◦ × × × × × × ×

✓ × × ✓ × ✓ ✓ ×

✓ ✓ ✓ ✓ ✓ × ✓ ✓

TABLE I. Comparison of quantum dot simulation frameworks (QDSim [18], QArray [15], QDFlow [10], QDarts [17], SimCATS [25], RF-Squad [16]). Here, “tunnel coupling” denotes coherent interdot coupling that hybridizes charge states and broadens CSD transition lines, as opposed to microscopic lead-barrier transport, and “latching” denotes metastable nonequilibrium charge retention on scan timescales. The symbol ◦ indicates a phenomenological implementation that reproduces the qualitative feature without the underlying physics: for latching, the simplified relaxation model originally introduced in QArray, which neglects the gate-voltage dependence of the relaxation rate and the possibility of relaxation through intermediate charge states; for tunnel coupling, the geometric rounding and interpolation of interdot transitions in SimCATS, which yields broadened transition lines and fractional occupations without a coherent Hamiltonian. See Appendix B for a comparison of QArray with the old latching implementation and QArray+.

-2.0 0.0 -2.0 (d) QArray+ Hubbard 376 ms 0.0

-2.0 -2.0

QArray+ Stochastic

41 ms

V1 (V) QArray+ Lindblad

0.0

V1 (V)

0.0

787 ms

V2 (V)

V1 (V)

V2 (V)

-2.0 -2.0 (c) 0.0

(b) 141 ms 0.0

V2 (V)

QArray

V2 (V)

(a) 0.0

V1 (V)

-2.0 0.0 -2.0

FIG. 1. Double dot charge stability diagram computed using four different models. a) The steady state charge stability diagram computed using QArray. b-d) The latched charge stability diagram computed using the QArray+ stochastic, spinless Hubbard and open quantum system models, respectively. The stochastic CSD shows latching but not hybridization of charge states. The spinless Hubbard model shows just hybridization of charge states. The open quantum system model shows both effects. The CPU compute time for each 200 × 200 CSD is indicated in the top-right corner.

plementary dynamical models for computing CSDs, each with differing trade-offs between computational complexity and the captured physical phenomena: - Stochastic Capacitance Model: This approach treats charge dynamics as classical Markov jump processes between discrete integer occupation states. Within this framework, the model effectively captures non-equilibrium latching phenomena, essential for simulating regimes where transition rates are slow relative to readout timescales. - Spinless Hubbard Model: To account for quantum mechanical effects in the steady state, this model implements a spinless Hubbard Hamiltonian. This allows for the computation of CSDs that exhibit coherent interdot hybridization, though it assumes a steady-state equilibrium that precludes the observation of non-equilibrium latching, mirroring the capabilities introduced in Ref. [16] (see Table I). - Open-Quantum-System Model: To provide a unified simulation of both coherent and incoherent processes, this model describes the array with a Lindblad master equation combining phonon and lead dissipators, which is evolved as a stochastic quantum-jump (Monte Carlo wavefunction) trajectory in the instantaneous eigenbasis. By bridging the gap between the previous two approaches, it enables the simultaneous treatment of coherent hybridization and the incoherent tunneling processes that lead to latching, offering a complete description of the transition between these regimes. Table I summarizes the capabilities of existing simulators along with those of QArray+, while Figure 1 shows

3 Vg1

Vg2 Cdg,11

n1

lead Γlead,1

tunnel

II. Cdg,22

n2

Cdd,12 t12 , Γinterdot,12 capacitive

MODELS FOR SIMULATING CHARGE STABILITY DIAGRAMS

lead

In this section, we describe the three dynamical models used to compute the charge stability diagrams in QArray+.

Γlead,2 A. tunnel and capacitive

FIG. 2. Schematic of a quantum dot device with tunnel and capacitive couplings. Each quantum dot carries charge ni and is tunnel-coupled to the source and drain leads with rates Γlead,i . The dots are mutually coupled through the inter-dot capacitance Cdd,12 = Cdd,21 , the coherent inter-dot tunnel coupling t12 , and the incoherent inter-dot transition rate Γinterdot,12 . Gate electrodes at voltages Vgi are capacitively coupled to the dots; for clarity, any cross-capacitances are omitted. The legend indicates the symbols used to represent tunnel, capacitive, and combined tunnel-capacitive couplings. Figure adapted from Fig. 1 of [26].

representative CSDs computed with the original QArray and with the new QArray+ models. In QArray+, all capacitance matrices and tunnel couplings include explicit gate-voltage dependence, allowing the simulator to reproduce device behavior across a wide range of operating conditions. These capabilities are provided through a unified interface, with optional functionality for generating charge-sensor signals and adding realistic noise to synthetic CSDs. The full stack is implemented in JAX, enabling just-in-time compilation and accelerator execution. In particular, the simulator is designed to scale from a single GPU to multi-GPU and multi-node environments. See Figure 2 for a schematic representation of a quantum dot device and Figure 3 for an overview of the structure of the software packages. In addition to these three dynamical models, the package provides a fourth simulation head—an integer groundstate (steady-state) solver that deterministically selects the minimum-energy charge configuration at each gatevoltage point, reproducing the steady-state model of the original QArray [used for the steady-state panel of Figure 1(a)]—alongside the original QArray model for backward compatibility as seen in Appendix B. In the following section II, we introduce the models used to compute CSDs. Then in section III we give examples of the CSDs generated by these models and discuss their features. Finally, in section IV we benchmark the models on CPUs, on a single GPU and on multiple GPUs.

Stochastic Capacitance Model

To model the non-equilibrium charge dynamics responsible for latching behavior, we evolve the charge configuration ⃗n of a quantum-dot array in discrete time steps, each corresponding to the integration time of the sensor measurement (i.e., one pixel of a charge-stability diagram). As a first approximation, we assume that within each integration window, at most one charge transition may occur, consistent with the experimentally relevant regime in which tunnel rates are low compared with the sampling bandwidth. As shown in section II C, a classical Monte Carlo model of exactly this structure emerges from the open-quantum-system treatment in the limit of vanishing interdot tunnel couplings, grounding the stochastic approach within the hierarchy of timescales that characterizes charge dynamics in quantum dot arrays. At the start of each time step, the algorithm enumerates all elementary charge-transfer processes: 1. Loading/Unloading: an electron tunnels from a reservoir into or out of a dot i, changing ni → ni ±1; 2. Interdot tunneling: an electron tunnels between neighboring dots i and j, conserving total electron number. For each possible transition, the change in electrostatic energy ∆U is computed according to the capacitance model for quantum dots [15–18, 26], such that U (⃗n) =

1 ⃗ T −1 ⃗ Q Cdd Q, 2

⃗ := |e|⃗n − Cdg V ⃗g , where Q

(1)

where matrices Cdd and Cdg denote the dot-dot and gate-dot capacitances, respectively, and may be functions of gate voltage. As in the constant-interaction picture, Equation 1 contains no quantized orbital level spacing. Although occupation-dependent orbital energies are not explicitly modeled, their smooth gate-voltage dependence can be partially emulated through the voltagedependent capacitance matrices [27]. That is, in general, ⃗g ) and Cdg := Cdg (V ⃗g ). Cdd := Cdd (V

(2)

In the following, we describe how tunneling rates are calculated for loading, unloading, and interdot transitions. a. Lead-mediated loading and unloading - Electron exchange between dot i and its reservoir occurs in the weak-tunneling regime and is therefore governed by Fermi–Dirac statistics. The loading rates are Γloading = i

4 ⃗g ), Cdg (V ⃗g ), Γ(V ⃗g ), tij (V ⃗g ) Gate-voltage dependence: Cdd (V Inputs ⃗g , device params V Cdd , Cdg (+ cross-caps) Electrostatic core (QArray) ⃗ = |e|⃗n − Cdg V ⃗g Q 1 ⃗ ⊤ −1 ⃗ U (⃗n) = 2 Q Cdd Q (gate-dependent optional)

(1) Integer ground state (steady-state) continuous ncont via convex QP; enumerate nearest; evaluate F (⃗n); JIT + VMAP batch (2) Stochastic latching Markov jumps per pixel; finite Γ, T , PSB (3) Tunnel-coupled ground state ⃗g ) + Ĥt (V ⃗g ); eigensolve Ĥ = ĤF (V

Outputs n (GS / latched / hybridized / Lindblad) Charge sensor Capacitive coupling; noise (1/f ); CSD map

(4) Lindblad open system P ρ̇ = −i[Ĥ, ρ] + e D[L̂e ]ρ; RK/MCWF Implementation: JAX + JIT → CPU / GPU / multi-GPU / multi-node

FIG. 3. Overview of QArray+ relative to QArray: shared electrostatic core with four interchangeable simulation heads (stochastic latching, tunnel-coupled ground state, Lindblad open-system evolution, and integer ground-state via continuous relaxation + candidate enumeration + energy selection), plus shared gate-voltage dependence, readout/noise, and JAX GPU acceleration.

⃗g )/(1 + exp(β∆Ui )), and the corresponding unloadΓi ( V ⃗g )/(1 + exp(−β∆Ui )), ing rates are Γunloading = Γi (V i where ∆Ui denotes the change in electrostatic energy associated with the ith loading transition, Γi is a (possibly gate-dependent) tunneling rate specified by the user, and β = (kB T )−1 is the inverse temperature. b. Interdot transitions - For interdot tunneling, the user specifies a rank-3 tensor Γinter s,ij , which parameterizes the transition rates between quantum dots i and j. The index s ∈ {0, 1} indicates whether Pauli spin blockade (PSB) is relevant for the transition: Γinter 0,ij corresponds to transitions that are not subject to PSB, whereas Γinter 1,ij corresponds to transitions that are subject to PSB. For a transition between dots i and j, the PSB rate is applied when both dots contain an odd number of charges prior to the transition. Similarly to the lead loading, we modulate the rate of interdot tunneling with the Fermi function to account for temperature effects, such that Γinter → Γinter /(1 + exp(β∆Uij )), where ∆Uij denotes ij ij the change in electrostatic energy upon an interdot transition between dots i and j, computed according to Equation 1. This Fermi factor is a bounded interpolation: it preserves the detailed-balance ratio e−β∆Uij between forward and backward transitions—and hence the correct equilibrium populations—while remaining finite at the interdot resonance, where the bosonic (phonon) factor of the open-system model (section II C) is instead resonantly enhanced, a divergence that is physically cut off by the phonon spectral density. This treatment is also an approximation in that it does not account for co-

herent tunneling between the dots; coherent tunneling is captured by the Lindblad-based open-system simulator discussed later. c. Stochastic charge state update - Having defined the transition rates, we now describe how the charge state is updated stochastically to reflect charge transfer events. For an integration window of duration τ , the probability that a transition with rate Γk occurs within this window is given by pk = 1 − exp(−Γk τ ).

(3)

To update the charge configuration, we perform a stochastic (dynamic Monte Carlo) update based on the transition probabilities defined above. The update procedure is as follows: 1. all allowed transitions are traversed in a random order; 2. for each transition, a Bernoulli trial with success probability pk is performed; 3. the first transition for which the trial succeeds is executed immediately; if no trial succeeds, the charge configuration remains unchanged. The procedure is then repeated for the next pixel of the charge stability diagram, starting from the updated charge configuration. This single-transition-per-time-step scheme is accurate in the slow-tunneling regime, where Γk τ ≪ 1 for all relevant transitions. The scheme reproduces the total jump

5 P Q probability exactly, since k (1 − pk ) = e−τ k Γk , but the random traversal order biases the selection among competing P channels relative to the exact branching ratios Γk / j Γj at second order in Γk τ [cf. Equation 32]. When multiple transitions within a single integration time must be resolved, or when Γk τ ≪ 1 is not satisfied, the method can be generalized by subdividing the integration window into Nr sub-intervals of duration τ /Nr and applying the update rule recursively within each subinterval, which suppresses this bias accordingly. d. Computational complexity - For a system consisting of ndot dots there are only Nr ndot (ndot + 1) charge transitions to evaluate, so the computational complexity of this approach scales polynomially with the number of dots. This polynomial scaling represents a dramatic improvement over steady-state solvers, which typically must evaluate the electrostatic energy of an exponentially growing number of charge states. See Appendix A for a complete discussion of the complexity.

B.

Tunnel coupling

The capacitance model used in the original QArray and in the previous model treats each charge configuration dot ⃗n ∈ Zn≥0 as a classical state with a well-defined electrostatic energy. However, in many experimentally relevant regimes, finite interdot tunnel coupling hybridizes charge configurations close to degeneracy, producing avoided crossings and fractional charge expectation values [28]. To address this limitation, we extend the capacitance model into a spinless Hubbard model, with a coherent tunneling Hamiltonian defined directly in the charge (occupation-number) basis, similar to that introduced in [16]. Then we compute CSDs by assuming that, for each gate voltage, the system adopts the lowest-energy eigenstate of this new charge-basis Hamiltonian. a. Charge-basis Hamiltonian - We work in the occupation-number basis {|⃗n⟩} = {|n1 , . . . , nndot ⟩} with ni ∈ {0, 1, . . . }. The diagonal (electrostatic) contribution is taken to be X ⃗g ) = ⃗g ) |⃗n⟩⟨⃗n| , ĤF (V F (⃗n; V (4) ⃗ n

⃗g ) is the electrostatic energy obtained from where F (⃗n; V the capacitance model, Equation 1. When the capacitances are themselves gate-voltage dependent, F should be read as a per-pixel electrostatic energy landscape rather than a thermodynamic free energy. In this formalism, overall constant prefactors (e.g. |e|2 /2) are absorbed into the choice of energy units; only relative energies enter the ground-state problem. Coherent interdot tunneling is modeled by single-mode hopping terms written in terms of charge ladder operators âj , â†i at each site act√ ni |. . . , ni − 1, . . .⟩ and ing as âi |. . . , ni , . . .⟩ = √ â†i |. . . , ni , . . .⟩ = ni + 1 |. . . , ni + 1, . . .⟩. With an

optional, user-specified gate-dependent tunnel matrix ⃗g ), the tunneling Hamiltonian is t(V X ⃗g ) = − ⃗g ) ↠âi . (5) Ĥt (V tij (V j i̸=j

Hermiticity of Ĥt requires tij = t∗ji ; in practice, the tunnel couplings are taken to be real and symmetric. In the charge basis, the only nonzero off-diagonal matrix elements connect states that differ by moving one electron from dot i to dot j: q ⃗g ) ni (nj + 1) if ⃗n′ = ⃗n − êi + êj , ⟨⃗n′ | Ĥt |⃗n⟩ = −tij (V (6) and vanish otherwise, where êi denotes the ith standard basis vector. Note that the charge ladder operators endow thephopping element with a soft-core bosonic enhancement ni (nj + 1) relative to the bare element −tij of a spinless fermionic Hubbard model with ni ∈ {0, 1} [16, 17]; the two conventions coincide for transitions between empty and singly occupied dots and differ only when multiply occupied dots participate. The full Hamiltonian is then ⃗g ) = ĤF (V ⃗g ) + Ĥt (V ⃗g ). Ĥ(V

(7)

b. Tunnel-coupled ground state and charge expecta⃗g , the tunneltion values - For each gate-voltage point V ⃗ coupled ground state |ψ0 (Vg )⟩ is obtained by solving the ⃗g ) |ψ0 ⟩ = E0 |ψ0 ⟩ and selecting the loweigenproblem Ĥ(V est eigenvalue. The simulator reports the expected dot occupations X ⃗g ) = ⟨ψ0 (V ⃗g )| n̂k |ψ0 (V ⃗g )⟩ = ⃗g )|2 nα,k , ⟨nk ⟩(V |ψ0,α (V α

(8) where α indexes the charge-basis states included in the calculation and nα,k denotes the occupation of dot k in basis state α. c. Efficient truncation of the charge basis - A naive implementation would require diagonalizing Ĥ in the full charge basis, whose dimension grows exponentially with ndot and the maximum charge per dot. Instead, QArray+ constructs a truncated charge basis of fixed size ntruncate that captures the low-energy manifold relevant near the ground state. ⃗g we first compute a continuConcretely, for each V ⃗g )—the minimizer of ous “open-charge” solution ⃗ncont (V Equation 1 with the integrality constraint on ⃗n relaxed— by solving the corresponding convex quadratic program. We then define an integer baseline ⃗nbase by taking the element-wise floor of ⃗ncont and generate a mixed-radix set of candidate integer configurations of the form ⃗n = ⃗nbase + ⃗δ, where δi is drawn from a small dot-dependent set. In the default heuristic used in the code, dots that couple most strongly to the gates are assigned a larger neighborhood δi ∈ {−1, 0, 1, 2}, more weakly coupled dots are assigned δi ∈ {0, 1}, and all remaining dots

6 are fixed by rounding; these asymmetric default neigh⃗g ) of borhoods are a heuristic choice. The energy F (⃗n; V each candidate is evaluated efficiently using a Cholesky factorization of Cdd , and the ntruncate lowest-energy candidates define the truncated basis. ⃗g ) is assembled in this truncated basis and Finally, Ĥ(V the ground state is computed either by a dense eigensolve (for small ntruncate ) or by a sparse Lanczos iteration (for larger truncated dimensions). This yields an accurate approximation to the tunnel-coupled ground state while avoiding the exponential scaling of the full basis. d. Computational complexity - Consider a system of ndot dots, and let ntruncate denote the truncated-basis dimension used for the eigensolve. For each gate-voltage ⃗g , the solver consists of: (i) generating a candipoint V ⃗g ) of integer charge configurations near the date set C(V ⃗g ), (ii) evaluating F (⃗n; V ⃗g ) for all relaxed solution ⃗ncont (V ⃗n ∈ C and selecting the ntruncate lowest-energy candidates, (iii) assembling the truncated Hamiltonian, and (iv) computing its ground state. ⃗g )| denote the number of candiLet Ncand := |C(V date configurations and Et the number of nonzero interdot tunnel couplings tij (i.e. edges in the tunnel graph). Given a Cholesky factorization of Cdd , each energy eval⃗g ) can be performed using two triangular uation F (⃗n; V solves and a dot product, scaling as O(n2dot ). Thus candidate scoring scales as O(Ncand n2dot ) per voltage point (up to an additional O(n3dot ) per point if Cdd must be refactorized due to gate-voltage dependence). The truncated Hamiltonian is sparse because Ĥt connects only charge states that differ by a single hop; the number of nonzero off-diagonal elements scales as nnz(Ĥt ) = O(ntruncate Et ). Computing the ground state by dense diagonalization costs O(n3truncate ), whereas a Krylov/Lanczos method costs O(niter nnz(Ĥ)) with niter iterations. Overall, the full scan therefore scales as    O N 2 Ncand n2dot + ntruncate Et + eig(ntruncate ) , where N is the linear dimension of the N × N gate-voltage grid and eig(ntruncate ) ∈ {O(n3truncate ), O(niter ntruncate Et )} depending on the eigensolver. Crucially, by keeping ntruncate (and in practice Ncand ) fixed and small, the computation is controlled by the truncated manifold rather than by the full charge basis, whose dimension grows exponentially with ndot . C.

Open quantum system

The stochastic charge-transition model described above treats charge dynamics as a classical Markov jump process between integer occupation states, capturing latching but not hybridization, while the spinless Hubbard model captures hybridization but not latching [29, 30]. To simulate both effects simultaneously, we

treat the quantum dot array as an open quantum system. To model the charge stability diagram while capturing both coherent hybridization and stochastic charging and relaxation dynamics, we employ an instantaneouseigenbasis stochastic Schrödinger equation (SSE) approach. Unlike standard Fock-basis quantum-jump models, which project the system onto unhybridized charge configurations at each jump, this formulation allows coherent superpositions to survive at avoided crossings between dissipative events [31]. a. Charge basis, Hamiltonian, and master equation - We again work in the charge basis {|⃗n⟩}, but now explicitly truncate the local occupations to ni ∈ {0, 1, . . . , nmax }, yielding a Hilbert-space dimension S = ⃗g , we con(nmax + 1)ndot . For each gate-voltage point V struct the same tunnel-coupled Hamiltonian ⃗g ) = ĤF (V ⃗g ) + Ĥt (V ⃗g ), Ĥ(V

(9)

with ĤF diagonal in the charge basis and Ĥt given by the hopping term in Equation 5. In this open-system setting, Ĥ generates coherent evolution between charge configurations that are connected by interdot tunneling. The system evolves under the Lindblad master equation (setting ℏ = 1) ⃗g ), ρ] + D(ρ), ρ̇ = −i [Ĥ(V

(10)

where the dissipator D comprises contributions from a bosonic phonon bath and from the fermionic leads: D(ρ) = Dphon (ρ) + Dlead (ρ).

(11)

b. Bosonic phonon bath - Incoherent interactions with phonons induce transitions between the instanta⃗g ), with energies Eα , while neous eigenstates |ψα ⟩ of Ĥ(V P preserving the total electron number N̂ = i n̂i . Because Ĥ conserves N̂ , every eigenstate carries a definite total charge Nα , and these processes drive relaxation and dephasing within each fixed-charge sector: X √  (12) Dphon (ρ) = D καβ |ψα ⟩⟨ψβ | (ρ), α̸=β

where D[L̂](ρ) = L̂ρL̂† − 21 {L̂† L̂, ρ} is the Lindblad dissipator and the sum runs over pairs of eigenstates within the same total-charge sector (Nα = Nβ ). The transition rates are those of a thermal bosonic bath, (   Γphon 1 + nB (|∆Eαβ |) , ∆Eαβ < 0, καβ = (13) Γphon nB (∆Eαβ ), ∆Eαβ > 0, where ∆Eαβ = Eα − Eβ , Γphon is the phonon coupling strength, and nB (E) = (eE/(kB T ) − 1)−1 is the Bose–Einstein occupation: downward transitions proceed at the stimulated-plus-spontaneous emission rate and upward transitions at the absorption rate. The rates obey the quantum detailed-balance relation for a thermal bath, καβ /κβα = e−∆Eαβ /(kB T ) , which follows from

7 [1 + nB (E)]/nB (E) = eE/(kB T ) . This bath acts as a dephasing mechanism, decohering superpositions on a timescale ∼ 1/Γphon while allowing relaxation within the fixed-charge manifold. c. Fermionic lead coupling - Charge carriers tunnel between the dots and external fermionic reservoirs, mediating transitions between charge states ⃗n and ⃗n ± êi , where êi is the ith standard basis vector. The corresponding dissipator is  X q D Dlead (ρ) = κi,± (⃗n) |⃗n ± êi ⟩⟨⃗n| (ρ), (14) i,±,⃗ n

where the sum retains only moves that respect the truncation 0 ≤ ni ≤ nmax . The transition rate is thermally regulated by the Fermi factor κi,± (⃗n) =

Γi  , 1 + exp ∆Ei,± (⃗n)/(kB T )

(15)

where Γi is the dot–lead coupling strength and ⃗g ) − U (⃗n; V ⃗g ) is the electrostatic ∆Ei,± (⃗n) = U (⃗n ± êi ; V energy cost of adding (+) or removing (−) an electron on dot i, computed according to Equation 1. These rates coincide with the lead-mediated loading and unloading rates used in the stochastic capacitance model. d. Hierarchy of timescales - Physical devices exhibit a clear separation of timescales, ∆0 ≫ Γphon ≫ Γlead ,

(16)

each eigenstate. For our charge-stability-diagram simulations, we initialize the first pixel in the instantaneous Hamiltonian’s ground state and then propagate this state using the stochastic wavefunction propagation described next. e. Stochastic wavefunction propagation - The dynamics of an open quantum system subject to inelastic transitions and decoherence are generally described by a master equation governing the density matrix ρ(t), but solving Equation 10 directly becomes computationally expensive for systems with many eigenstates. An alternative approach is to represent the evolution stochastically using the quantum-jump method (also known as Monte Carlo wavefunction simulation): pure-state trajectories evolve under a non-Hermitian effective Hamiltonian, interspersed with stochastic collapse events representing inelastic transitions. When averaged over many trajectories, this approach reproduces the master equation’s predictions while remaining numerically efficient. In our case, however, we are not interested in the ensemble average, as charge stability diagrams are typically measured in a small number of shots; a single trajectory directly mimics an individual experimental scan. Under the quantum-jump protocol, the state |Ψ(t)⟩ evolves over a timestep ∆t = τint according to the following procedure. We first expand the state in the in⃗g ): stantaneous eigenbasis of Ĥ(V X |Ψ(t)⟩ = cβ (t) |ψβ (t)⟩ . (18) β

where ∆0 ∼ max |spec(Ĥ)| sets the overall spectral scale of Ĥ and Γlead denotes the typical magnitude of the dot– lead couplings Γi . This ordering partitions the dynamics into three distinct regimes. On the coherent timescale τcoh ∼ 1/∆0 , the system exhibits nearly coherent evolution between electronic eigenstates, with dissipation negligible. On the phonon timescale τphon ∼ 1/Γphon , the phonon bath gradually decoheres oscillations through incoherent transitions that preserve charge. Finally, on the lead timescale τlead ∼ 1/Γlead , charge transfer between dots and reservoirs modulates the total charge state. When measuring a charge stability diagram, we assume that the gate voltages change instantaneously between adjacent pixels and that the charge-sensor signal is integrated over a time interval τint , corresponding to one pixel. We simulate the experimentally realistic regime τphon ≪ τint ≲ τlead ,

(17)

ensuring that the integration time is much longer than phonon relaxation timescales but comparable to the leadinduced decoherence timescale; for reference, typical integration times are on the microsecond timescale. Under these conditions, coherences between eigenstates are damped out within each pixel, and we treat the system as evolving stochastically among the Hamiltonian’s instantaneous eigenstates. Note that because these eigenstates are themselves hybridized superpositions of charge configurations, interdot charge coherence is retained within

The non-Hermitian effective Hamiltonian iX Ĥeff = Ĥ − Λα |ψα ⟩⟨ψα | 2 α

(19)

generates the unnormalized coefficients over the timestep ∆t:    1 cβ (t + ∆t) = cβ (t) exp −iEβ − Λβ ∆t , (20) 2 P where Λβ = α̸=β Rαβ is the total escape rate from eigenstate |ψβ ⟩, with Rαβ a matrix encoding the transition rate from eigenstate |ψβ ⟩ to eigenstate |ψα ⟩; we outline how the matrix Rαβ is populated below. The imaginary part of the exponent encodes the coherent energy evolution, while the real part accounts for exponential decay due to dissipation—this is the non-Hermitian evolution that underlies the master-equation approach. Because the gate voltages are constant within a pixel, the eigenbasis is fixed during each timestep, and Equation 20 follows exactly from Equation 19. The probability that at least one quantum jump occurs during the timestep is X Pjump = 1 − |cβ (t + ∆t)|2 . (21) β

At each step, a uniform random variable ξ ∈ [0, 1) determines which of two scenarios occurs.

8 Case 1: No jump (ξ ≥ Pjump ). The state undergoes continuous evolution under Ĥeff , and the coefficients are renormalized to restore unit norm: P β cβ (t + ∆t) |ψβ (t + ∆t)⟩ p . (22) |Ψ(t + ∆t)⟩ = 1 − Pjump Case 2: Quantum jump (ξ < Pjump ). A stochastic collapse occurs, representing an inelastic transition. The source eigenstate |ψj ⟩ is sampled according to the probability that the norm loss originated from that state,
|cj (t)|2 1 − e−Λj ∆t P (source = j) = , (23) Pjump which reduces to P (source = j) ∝ |cj (t)|2 Λj to first order in ∆t. Given the source, the target eigenstate |ψi ⟩ is randomly selected via the branching ratio P (target = i | source = j) =

Rij . Λj

(24)

The state then collapses to the target eigenstate: |Ψ(t + ∆t)⟩ = |ψi (t + ∆t)⟩ .

(25)

Note that, mirroring the single-transition-per-time-step scheme of the stochastic capacitance model, at most one quantum jump is applied per evolution substep (by default, a single substep spans the whole integration window). In parallel mode (introduced below), each row is re-initialized in the instantaneous ground state, so this one-jump update acts on a state that has already relaxed. In continuous mode, where the trajectory is carried between pixels, the regime of Equation 17 implies Γphon τint ≫ 1, and the integration window should therefore be subdivided into Nr ≳ Γphon τint substeps—in analogy with the recursion depth of the stochastic capacitance model—so that the jump probability per substep remains small. The simulated charge signal for a given pixel is then obtained by computing the expectation value of the charge ˆ |Ψ(t + ∆t)⟩, where the occupation vector, ⟨Ψ(t + ∆t)| ⃗n ˆ components of ⃗n = (n̂1 , . . . , n̂ndot ) correspond to the number operators of the individual quantum dots. Because the eigenstates are hybridized superpositions of charge configurations, this signal captures the fractional charge expectation values responsible for avoided crossings, while the stochastic jumps reproduce latching dynamics; in the limit of vanishing tij , the diagrams reduce to the classical latching picture. f. Scan modes - When computing a twodimensional charge stability diagram, the gate voltages are typically swept in a raster pattern. To simulate this, we provide two distinct computational modes. In continuous mode, the entire diagram evolves sequentially: the final charge configuration of one pixel serves as the initial state for the subsequent pixel, preserving the stochastic history of the system across the entire

scan, including the flyback between rows. Alternatively, in parallel mode, we assume the system fully relaxes during the raster flyback or the subsequent wait time before starting the row, such that each new row is initialized in the instantaneous ground state. This assumption removes the physical dependence between rows, rendering the calculation embarrassingly parallel and significantly reducing computation time. g. Dissipative jump rates in the eigenbasis - To simulate the stochastic trajectories described above, we must compute the transition-rate matrix Rαβ , which defines the probability per unit time of a quantum jump from the instantaneous source eigenstate |ψβ ⟩ to the target eigenstate |ψα ⟩. The total transition rate is the sum of the independent contributions from the phonon bath and the fermionic leads: phon lead Rαβ = Rαβ + Rαβ .

(26)

Crucially, by constructing the jump rates directly between the instantaneous energy eigenstates, we are applying the secular approximation, thereby explicitly dropping the jump-induced coherences (off-diagonal terms in the density matrix). This approximation is physically justified by the system’s separation of timescales (∆0 ≫ Γphon , Γlead ), ensuring that the coherent energy scale strictly dominates and any rapidly oscillating offdiagonal coherences average to zero. As a result, inelastic collapse events map populations directly to populations in the instantaneous eigenbasis. We note that ∆0 characterizes the overall spectral scale; the locally relevant frequency for a given pair of eigenstates is their splitting |Eα −Eβ |, and the secular condition therefore fails locally at avoided crossings, where the gap collapses to ∼ 2tij and can become comparable to the dissipative rates. In these narrow regions the secular construction should be regarded as an approximation. The phonon dissipator Dphon is already defined via transitions between system eigenstates. Thus, the phonon transition rate follows directly from the detailedbalance coefficients of Equation 13, phon Rαβ = καβ

for Nα = Nβ and α ̸= β,

(27)

where Nα and Nβ are the total charges of the target and source eigenstates, respectively; all other elements phon vanish, and in particular Rαα = 0. For the lead dissipator Dlead , the jump operators are defined in the local charge basis |⃗n⟩. To find the eigenstate-to-eigenstate transition rate, we must project these jump operators into the eigenbasis. Let an instantaneous eigenstate be expanded in the charge baP (γ) (γ) sis as |ψγ ⟩ = n⟩, where c⃗n = ⟨⃗n|ψγ ⟩. Each ⃗ n c⃗ n |⃗ Lindblad jump operator L̂ contributes a transition rate | ⟨ψα | L̂ |ψβ ⟩ |2 from |ψβ ⟩ to |ψα ⟩. Summing the contri-

9 butions of all single-electron lead jump operators yields q  X 2 lead ⟨ψα | κi,± (⃗n) |⃗n ± êi ⟩⟨⃗n| |ψβ ⟩ Rαβ =

In the absence of coherent evolution, the probability that at least one charge-transfer event occurs during the integration window is simply

i,±,⃗ n

=

X

(α)

κi,± (⃗n) c⃗n±êi

2

(β) 2

c⃗n

.

(28)

i,±,⃗ n

This result carries a clear physical interpretation: the rate of an electron tunneling into or out of the dot array is governed by the underlying tunneling rate κi,± (⃗n), (β) weighted by the probability |c⃗n |2 that the initial eigenstate is in the necessary pre-tunneling charge configu(α) ration |⃗n⟩, and multiplied by the probability |c⃗n±êi |2 that the target eigenstate supports the resulting posttunneling charge configuration |⃗n ± êi ⟩. Summing the (i, ±, ⃗n) channels incoherently drops interference between different charge-basis pathways relative to a single coherent lead operator; this is the natural choice here because the rate κi,± (⃗n) depends explicitly on the pre-jump configuration ⃗n. Finally, the total escape rate from a given eigenstate |ψβ ⟩ is obtained by summing over all possible target states α: X Λβ = Rαβ . (29) α̸=β

This rate directly dictates the non-Hermitian decay and the total jump probability for the stochastic wavefunction propagation. h. Classical limit: stochastic capacitance model - In the limit where interdot tunneling vanishes (tij = 0 for all i, j), the eigenstates correspond exactly to the charge Fock states |⃗n⟩. This allows for considerable computational simplification: the state of the system is fully captured by the discrete charge vector ⃗n, and the stochastic wavefunction propagation reduces to a classical Monte Carlo simulation over the charge space. In this limit, for a given charge state ⃗n, the fermionic jump operators defined in Equation 14 permit transitions to the 2ndot charge states differing by one charge in a single dot, with rates computed according to Equation 15. Similarly, the bosonic jump operators permit internal charge rearrangements that conserve the total electron number. Because the number of such charge-preserving configurations scales exponentially with the number of dots, we restrict our attention to the ndot (ndot − 1) pairwise interdot transitions, in which a single electron is swapped between two dots while the charge state of all other dots remains unchanged; the corresponding phonon-induced transition rates are computed according to Equation 13. To evolve the system over an integration window of duration τint , we apply a stochastic jump protocol that directly mirrors the wavefunction jump procedure. First, we compute the total escape rate from the current charge state ⃗n by summing over all allowed fermionic and bosonic transition rates Rk : X Λ⃗n = Rk . (30) k

Pjump = 1 − exp(−Λ⃗n τint ).

(31)

A uniform random variable ξ ∈ [0, 1) determines the outcome for the given pixel. If ξ ≥ Pjump , no jump occurs: the system remains in the current state and the charge configuration ⃗n is unchanged. If ξ < Pjump , a charge-transfer event occurs, and the specific transition k is randomly selected from the set of allowed transitions according to its branching ratio: P (target = k | jump) =

Rk . Λ⃗n

(32)

The charge configuration is immediately updated to the charge state associated with transition k. This procedure is then repeated for the next timestep of the charge stability diagram, starting from the resulting charge configuration. As with the stochastic wavefunction propagation, we implement continuous and parallel modes to compute the scan. While this method sacrifices the ability to capture the coherent hybridization of charge states, it is considerably more computationally efficient, scaling polynomially rather than exponentially with the number of dots. This limit recovers the structure of the stochastic capacitance model introduced in section II A, with the interdot rates now inherited from the phonon bath rather than specified directly by the user. The two interdotrate conventions share the same detailed-balance ratio, and hence the same equilibrium populations; they differ only in the lineshape of thermally broadened interdot transitions, where the Fermi factor of section II A saturates while the bosonic rate of Equation 13 is resonantly enhanced. i. Steady-state limit: spinless Hubbard model - We also consider the complementary regime of long integration times, in which the integration window exceeds even the lead timescale, such that τphon ≪ τlead ≪ τint .

(33)

Under these conditions, the quantum dot array fully relaxes at each gate-voltage configuration and, for temperatures small compared to the relevant energy gaps, settles into its instantaneous ground state. Consequently, every pixel of the charge stability diagram is physically independent and can be computed in parallel. We emphasize that Equation 17 and Equation 33 describe two distinct operating regimes of the same model—fast scans that exhibit latching versus slow scans that equilibrate— selected by the experimental integration time, rather than competing assumptions. In this limit, the opensystem model reduces to the tunnel-coupled ground-state (spinless Hubbard) treatment described above, including the truncated-basis construction that avoids diagonalizing the Hamiltonian in the exponentially large full charge basis.

10 j. Computational complexity - For the general opensystem evolution, the dominant cost is governed by the Hilbert-space dimension S = (nmax + 1)ndot : each gate⃗g ) voltage point requires the eigendecomposition of Ĥ(V and the assembly of the rate matrix Rαβ , both of which scale polynomially in S. The classical limit above avoids this cost entirely, which is why the stochastic capacitance model remains tractable for large arrays.

III. A.

RESULTS

Stochastic Capacitance Model

Figure 4 illustrates a series of hole-based double-dot charge stability diagrams computed using both the original QArray algorithm and the newly developed QArray+ approach. Each diagram is evaluated on a 200×200 gatevoltage grid, with the corresponding compute time shown in the top-right corner of each panel and the extracted charge state overlaid. Panel 4(a) shows the steadystate solution obtained using the original QArray method, which relies on numerically solving for the stationary occupation probabilities at each point in gate space. In contrast, panel 4(b) displays the charge stability diagram generated by the QArray+ stochastic solver in the idealized limit of infinite tunnel rates and zero lead temperature, where the system deterministically relaxes to the minimum-energy charge configuration. Panels 4(c)–(g) demonstrate latched charge stability diagrams under finite tunnel rates. In this regime, we assume that each horizontal sweep (each “row” of the voltage grid) begins after sufficient relaxation time such that the system returns to its ground-state charge configuration. This assumption enables efficient parallel computation across rows. Panel 4(h) shows the fully isolated limit, in which the dots are disconnected from their reservoirs. Here, the user specifies an initial charge state at the lower-left corner, and the state is stochastically propagated across the grid. As no relaxation occurs, each row depends on the evolution of the previous one, and parallelization is therefore not possible. All simulations in Fig. 4 use Nr = 1 and simulate the effect of a charge sensor incorporating 1/f noise. To explore the impact of increasing experimental realism, Fig. 5 presents simulated charge stability diagrams incorporating Pauli spin blockade (PSB) and finite-temperature effects. Panels 5(a–b) display a double-dot diagram computed at a finite temperature of kB T = 0.04, expressed in units of the average dot charging energy. In this simulation, the first dot is strongly coupled to the lead (Γ1 τ = 100), while the second dot is more weakly coupled (Γ2 τ = 1). Near the charge transition, the high tunneling rate to the lead causes the charge state to fluctuate rapidly between (1, 1) and (0, 1). While these individual stochastic events are resolved in the raw simulation, they are typically not observed experimentally due to the finite bandwidth of

the measurement setup; instead, the measured signal represents a time average over multiple transitions. In Panel 5(b), we simulate this experimental limit by increasing the recursion depth Nr to allow for multiple stochastic jumps within a single timestep. The resulting time-averaged signal yields a characteristic “doublestep” feature frequently observed in experimental measurements (see Appendix D). Panel 5(c) illustrates the effect of Pauli spin blockade, which suppresses the (1, 1) → (0, 2) transition and manifests as latching at the interdot transition boundary. Finally, Panel 5(d) demonstrates the model’s capability to incorporate gate-voltage dependence within the capacitance matrices and tunnel rates, further aligning the simulation with the non-linearities of physical device behavior.

B.

Spinless Hubbard model

The stochastic capacitance model captures nonequilibrium latching, but by construction treats each integer charge configuration as a classical state and therefore cannot reproduce coherent interdot hybridization. We next quantify the impact of coherent interdot tunneling in the equilibrium limit using the spinless Hubbard model described above. Figure 6 shows hole-based double-dot charge-stability diagrams computed at zero temperature for increasing interdot tunnel coupling t12 . In the limit t12 = 0 [Fig. 6(a)], the result reduces to the classical ground-state picture: the diagram is piecewise constant with sharp boundaries between integer charge states, yielding the familiar honeycomb pattern. For finite t12 [Fig. 6(b–d)], the interdot charge degeneracies are hybridized. Concretely, near a degeneracy between two configurations connected by a single interdot hop, e.g. |n1 , n2 ⟩ and |n1 − 1, n2 + 1⟩, the effective two-level sector exhibits an avoided crossing with gap ∼ 2t12 , and the ground-state occupations become fractional : ⟨n1 ⟩ and ⟨n2 ⟩ interpolate continuously across the interdot transition instead of switching discontinuously [27]. In the simulated sensor response this appears as a progressive broadening and rounding of the interdot transition lines and a smoothing of features around the triple-point regions as t12 increases. An important structural point is that the coherent tunnel Hamiltonian Ĥt conserves total charge. As a result, only charge configurations within a fixed totaloccupation manifold hybridize, so lead-transition boundaries (which separate different total-charge sectors) remain sharp in this equilibrium solver. This cleanly isolates the physics of hybridization, producing avoided crossings and charge delocalization free of any competing dissipative effects. It also motivates the open-system treatment below, which becomes necessary when finiterate exchange with reservoirs and scan-history effects are essential.

11

V2 (V)

(a) 0.5

steady state

(b) 118 ms

(2; 0)

(1; 0)

(0; 0)

(2; 1)

(1; 1)

(0; 1)

(2; 2)

(1; 2)

(0; 2)

-2.0 (e) 0.5

(1; 0)

(0; 0)

(2; 1)

(1; 1)

(0; 1)

(2; 2)

(1; 2)

(0; 2)

42 ms (2; 0)

(1; 0)

(0; 0)

(2; 0)

(1; 0)

(0; 0)

(2; 1)

(1; 1)

(0; 1)

(2; 2)

(1; 2)

(0; 2)

42 ms (2; 0)

(1; 0)

(0; 0)

(2; 1)

(1; 1)

(0; 1)

(2; 2)

(1; 2)

(0; 2)

(h)

~¡ ¿ = (1; 0)

42 ms (1; 0)

~¡ ¿ = (0:1; 1)

42 ms (2; 0)

(g)

~¡ ¿ = (0; 1)

(d)

~¡ ¿ = (1; 1)

41 ms

(2; 0)

(f)

~¡ ¿ = (1; 0:1)

(c)

~¡ ¿ = (1; 1)

~¡ ¿ = (0; 0)

42 ms (2; 0)

(1; 0)

(0; 0)

(2; 1)

(1; 1)

(0; 1)

(2; 2)

(1; 2)

(0; 2)

88 ms (3; 1)

V2 (V)

(0; 0)

(2; 1)

(1; 1)

(0; 1)

(2; 2)

(1; 2)

(0; 2)

-2.0 -2.0

V1 (V)

(2; 1)

(2; 2)

0.5 -2.0

(1; 1)

(1; 2)

(0; 1) (0; 2)

V1 (V)

0.5 -2.0

V1 (V)

(2; 2)

(1; 3)

0.5 -2.0

V1 (V)

0.5

FIG. 4. Hole-based double-dot charge stability diagrams, each computed on a 200 × 200 gate-voltage grid, with the compute time indicated in the top-right corner and the charge state overlaid. (a) Steady-state solution obtained using the original QArray algorithm. (b) Charge stability diagram computed with the QArray+ stochastic solver in the limit of infinite tunnel rates and zero lead temperature, such that the system deterministically adopts the lowest-energy charge configuration. (c–g) Latched charge stability diagrams for finite tunnel rates. It is assumed that sufficient time elapses between the end of one row and the start of the next for the system to relax to the ground charge state at the beginning of each row, allowing each row to be computed in parallel. (h) Charge stability diagram for the fully isolated case, where the dots are decoupled from the charge reservoirs, so that only interdot transitions are present. The user specifies the initial charge state in the lower-left corner, and this state is propagated stochastically across the diagram; in this regime, parallelization is not possible. All simulations use Nr = 1 and include a charge-sensor readout with 1/f noise.

C.

Open quantum system

Finally, we combine coherent hybridization with finiterate dissipative tunneling using the open-quantumsystem model. In this setting, the relevant competition is between (i) the coherent interdot hybridization scale set by t12 , (ii) the dissipative lead tunneling rates Γlead,i , and (iii) the effective dwell time per pixel τint . When the characteristic relaxation rate Γrelax set by these dissipative couplings is fast on the scale of a pixel (Γrelax τint ≫ 1), the diagram approaches the equilibrium (Hubbard) limit; when relaxation is slow (Γrelax τint ≲ 1), scan-history dependence emerges and latching-like metastability becomes visible even in the presence of coherent coupling. Figure 7 demonstrates this crossover using the Lindblad solver in a regime where the device is asymmetrically coupled to the reservoirs. We set the temperature to kB T = 0.01 (in units of the mean charging energy) and evolve the system for a fixed interval τint per pixel, propagating the quantum state along the scan direction in parallel mode, i.e., with each row initialized in the instantaneous ground state. To isolate the interplay between coherent hybridization and reservoir exchange,

phonon-mediated incoherent interdot transitions are disabled in this example (Γphon = 0), so redistribution between dots occurs only through the coherent tunnel coupling t12 . With the phonon bath switched off, the hierarchy of Equation 16 reduces to the single requirement ∆0 ≫ Γlead , which is satisfied here: the eigenbasis construction only requires the coherent scale to dominate all active dissipative rates. Panels 7(a,c) show the case where only dot 1 is coupled to a reservoir, Γlead = (1, 0), while panels 7(b,d) couple only dot 2, Γlead = (0, 1). When the lead is attached to dot 1, the system can efficiently load and unload charge through that dot; combined with coherent hybridization, this provides an effective relaxation pathway for the two-dot manifold and yields comparatively smooth, near-equilibrated transition features. In contrast, when the lead is attached to dot 2, charge exchange occurs through a different physical bottleneck relative to the scan axes, and the system can become trapped in metastable charge configurations over extended scan ranges, producing jagged, scan-direction-dependent transition features characteristic of latching. Increasing the coherent coupling from t12 = 0.05 to t12 = 0.10 [Fig. 7(b)

12

(0; 0)

(2; 1)

(1; 1)

(0; 1)

(2; 2)

(1; 2)

(0; 2)

(1; 0)

(0; 0)

(2; 1)

(1; 1)

(0; 1)

(2; 2)

(1; 2)

(0; 2)

(2; 0)

(2; 0)

(1; 0)

(0; 0)

(2; 1)

(1; 1)

(0; 1)

(2; 2)

(1; 2)

(0; 2)

-1.8 -1.8

V1 (V)

0.2 V1 (V) voltage dependent 84 ms

V2 (V)

V2 (V)

-1.8 -1.8 -1.8 0.2 -1.8 V1 (V) (c) ~¡ ¿ = (100; 1) and PSB (d) 0.2 155 ms 0.5

-4.0 0.2 -4.0

t12 = 0.00

-3.0 -3.0 (c) 1.0

V1 (V)

0.5

FIG. 5. Hole-based charge stability diagrams illustrating how increasing levels of experimental realism—including charge-sensor readout, Pauli spin blockade (PSB), and finite-temperature effects— modify the observed features. All simulations are performed on a 200 × 200 voltage grid. (a–b) Double-dot charge stability diagrams simulated at kB T = 0.04 for Nr = 1 and Nr = 10, respectively. Fast tunneling between the first dot and the leads, together with finite-temperature effects near the lead transitions, make loading and unloading comparably likely and lead to stochastic charge-state switching. For Nr = 1, this leads to unphysical artifacts in the CSD, while for Nr = 10 the fluctuations are captured correctly, resulting in the experimentally observed doubling of the transition. (c) Double-dot charge stability diagram simulated at kB T = 0.01, demonstrating the effect of Pauli spin blockade, which suppresses the (1, 1) → (0, 2) transition. (d) Double-dot charge stability diagram incorporating gate-voltage-dependent interdot capacitive coupling and lead tunnel rates.

vs. (d)] suppresses these non-equilibrium artifacts by increasing the hybridization-induced connectivity between relevant charge configurations, thereby accelerating relaxation through the lead-coupled dot and pushing the diagram toward the equilibrium limit. Because the Lindblad solver propagates the quantum state across the scan and samples stochastic quantum jumps, it is computationally more expensive than the equilibrium ground-state solver for the same grid resolution; nevertheless, the JAX implementation keeps these simulations tractable and enables accelerator execution for large-scale dataset generation.

-3.0 -3.0

(b) 570 ms 1.0

t12 = 0.01

350 ms

V2 (V)

(1; 0)

(a) 1.0

-3.0 1.0 -3.0 (d) t12 = 0.10 356 ms 1.0

V1 (V)

V1 (V) t12 = 0.20

1.0

355 ms

V2 (V)

(2; 0)

~¡ ¿ = (100; 1) 378 ms

V2 (V)

(b) 44 ms 0.2

V2 (V)

~¡ ¿ = (100; 1)

V2 (V)

V2 (V)

(a) 0.2

V1 (V)

-3.0 1.0 -3.0

V1 (V)

1.0

FIG. 6. Tunnel-coupled ground-state (spinless Hubbard) charge stability diagrams for increasing interdot tunnel coupling t12 . Coherent tunneling hybridizes near-degenerate interdot charge configurations, producing avoided crossings and fractional charge expectation values that broaden the interdot transition lines as t12 increases. Each panel shows a 200 × 200 gate-voltage scan; compute times on CPU are indicated in the top-right corners. The larger compute time for the t12 = 0 panel reflects one-time just-in-time compilation included in its timing; the remaining panels reuse the compiled kernel.

IV.

BENCHMARKING

All the CSDs displayed so far were computed using the CPU of an M3 MacBook Pro; the compute time on this hardware is superimposed in the top right-hand corner. We now benchmark the scaling of QArray+ to larger arrays and quantify accelerator speedups on a server-class system. All benchmarks were run on a node equipped with 8× NVIDIA H100 SXM5 80 GB GPUs and Intel Xeon Platinum 8480C CPUs. Unless stated otherwise, each point corresponds to the profiler time required to compute a full 100 × 100 charge-stability diagram (CSD) after JIT compilation and warm-up. Since each benchmark computes a fixed-resolution scan, the reported times scale approximately linearly with the number of pixels, allowing extrapolation to other grid sizes (e.g., 200 × 200 is ≈ 4× more work than 100 × 100) as seen in the Appendix. Figure 8 summarizes three representative workloads corresponding to the main simulation heads used in QArray+. We highlight the following subgraphs:

13

(b) 800 ms 0.5

t12 = 0.05 ~Γ τ = (0, 1)

790 ms

V2 (V)

t12 = 0.05 ~Γ τ = (1, 0)

V2 (V)

(a) 0.5

-2.0 -2.0

V2 (V)

V2 (V)

-2.0 -2.0 -2.0 0.5 -2.0 0.5 V1 (V) V1 (V) (c) t12 = 0.10 ~Γ τ = (1, 0) (d) t12 = 0.10 ~Γ τ = (0, 1) 0.5 785 ms 0.5 786 ms

V1 (V)

-2.0 0.5 -2.0

V1 (V)

0.5

FIG. 7. Open-system (Lindblad) charge stability diagrams combining coherent hybridization and finiterate tunneling to leads. The panels compare two interdot tunnel couplings (t12 = 0.05 and t12 = 0.10) and two asymmetric reservoir couplings, Γlead = (1, 0) and Γlead = (0, 1), illustrating the emergence and suppression of scan-historydependent latching as the balance between coherent tunneling, dissipative tunneling, and pixel dwell time is varied. Each panel shows a 200 × 200 scan; compute times on CPU are indicated.

a. Stochastic latching (capacitance Markov-jump solver). Fig. 8(a) compares the original QArray baseline implementation to the QArray+ latching solver on CPU, single GPU, and multi-GPU execution. The baseline shows a rapid growth in runtime with dot count and becomes impractical beyond ∼16 dots over the tested range, consistent with the exponential growth of the underlying state space in steady-state formulations. In contrast, the QArray+ latching solver exhibits the expected polynomial scaling in ndot and remains tractable out to 64 dots. At the largest system size shown, the CPU backend requires ∼2.8 s per 100×100 scan, while a single H100 reduces this to ∼0.4 s and 8×H100 reduces it further to ∼0.17 s, demonstrating substantial end-to-end speedups from accelerator execution. b. Tunnel-coupled ground state (spinless Hubbard head). Fig. 8(b) benchmarks the tunnel-coupled ground-state solver against QDarts. While QDarts becomes expensive already for ≳ 5 dots (hundreds of seconds per scan in this configuration) and extrapolates to prohibitive runtimes at larger ndot , QArray+ remains in the few-second to few-tens-of-seconds regime up to 20 dots. This improvement is enabled by the fixed-size truncated charge basis and the use of sparse linear

algebra (Lanczos-type eigensolvers), which avoid explicit enumeration and diagonalization in the full exponentially large charge basis. Across the tested range, sparse multi-GPU execution provides an additional speedup over a single GPU, with the benefit depending on the effective sparsity of the tunnel graph and the resulting per-point Hamiltonian construction and eigensolve cost. c. Open-system Markov jumps (Lindblad MCWF eigenjump head). Fig. 8(c) reports performance of the open-system Lindblad solver. The accelerator execution substantially extends the accessible regime of this fully physics-informed solver: at 10 dots, the CPU runtime reaches ∼1.3×103 s (∼22 min) per scan, whereas a single H100 reduces this to ∼110 s (∼2 min), and 8×H100 reduces it to ∼30 s. Overall, these benchmarks show that QArray+ supports high-throughput dataset generation across multiple physical regimes, and that JAX-based accelerator execution, especially multi-GPU sharding, yields orderof-magnitude reductions in end-to-end runtime for the workloads most relevant to large-scale simulation of quantum-dot arrays.

V.

CONCLUSION

We have introduced QArray+, a versatile simulator for modeling non-equilibrium charge dynamics in tunnelcoupled semiconductor quantum-dot arrays. By providing a quantum open-system description based on Lindblad master-equation evolution, the framework enables a unified treatment of coherent dynamics and dissipative tunneling to leads. This approach allows for the generation of charge-stability diagrams (CSDs) that simultaneously capture complex phenomena such as latching and interdot tunnel coupling. In regimes where coherent effects are less dominant, QArray+ allows for a transition to a stochastic model characterized by polynomial rather than exponential scaling. This efficiency enables the simulation of largescale quantum-dot arrays that would be computationally prohibitive using conventional methods. Consequently, we anticipate that this framework will serve as a robust “digital twin” for the offline testing and verification of large-scale devices. The framework naturally incorporates gate-voltagedependent capacitances, Pauli spin blockade, finitetemperature effects, realistic sensor readout and noise models. These features make QArray+ suited for integration into automated tuning pipelines and the generation of large-scale synthetic datasets for machine-learning applications. As quantum-dot devices grow in complexity and fast measurements become widespread, the ability to bridge tunnel-coupling-induced hybridization with dissipative non-equilibrium dynamics will be essential. The methods presented here lay the groundwork for scalable, realistic simulation tools designed to support the next generation of semiconductor quantum processors. Future

14 (a)

(b)

QArray vs QArray+ Stochastic

(c)

QArray+ Hubbard

QArray+ Lindblad

10 3

Compute time (s)

Compute time (s)

10 0

10 −1

10 −2

QArray CPU GPU multi-GPU

0

8

16

24

32

40

Number of dots

48

56

64

Compute time (s)

10 2

10 1

10 0

GPU multi-GPU QDarts

2

6

10

14

Number of dots

18

10 2

10 1

10 0

CPU GPU multi-GPU

2

3

4

5

6

7

Number of dots

8

9

10

FIG. 8. Benchmarks for QArray+ on latching, tunnel coupling and Lindblad features. All benchmarks use a 100 × 100 charge-stability-diagram resolution on an 8 × NVIDIA H100 (80 GB HBM3) node. (a) QArray baseline versus QArray+ latching on CPU, single GPU, and 8-GPU, showing that QArray+ on GPU and multi-GPU outperforms the baseline at all dot counts scaling up to 64 dots; multi-GPU maintains sub-second runtimes throughout. (b) Tunnel-coupling benchmark comparing QArray+ (single GPU and multi-GPU) against QDarts, which exceeds 200s at 5 dots while QArray+ remains tractable to 20 dots. (c) Lindbladian benchmark; multi-GPU provides an order-of-magnitude speedup over CPU at 10 dots.

work using QArray+ will also enable rapid data generation for training autonomous tuning pipelines, which are crucial for realizing scalable semiconductor quantum devices. ACKNOWLEDGMENTS

P.V. performed this research as part of an internship at NVIDIA. N.A. acknowledges support from the European Research Council (grant agreement 948932), the Royal Society (URF-R1-191150) and the United States Army Research Office under Award No. W911NF-24-20043. This research was also supported by the European Union through the Horizon 2020 research and innovation program under the Grant Agreement No. 951852 (QLSI) and the Horizon Europe Framework Program under grant agreement No. 101069515 (IGNITE). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union, Research Executive Agency or UK Research & Innovation. Neither the European Union nor UK Research & Innovation can be held responsible for them. P.V. is supported by the United States Army Research Office under Award No. W911NF-21-S-0009-2. A.P., B.v.S., and M.V. acknowledge support from the NWO through a Vici grant (VI.C.242.031) and the National Growth Fund program Quantum Delta NL (grant NGF.1582.22.001). This research was sponsored in part by the Army Research Office (ARO) under Award No. W911NF-23-1-0110 and by The Netherlands Ministry of Defense under Award No. QuBits R23/009. The views, conclusions, and recommendations contained in this document are those of the authors and are not necessarily endorsed nor should they be interpreted as representing the

official policies, either expressed or implied, of the Army Research Office (ARO) or the U.S. Government, or The Netherlands Ministry of Defense. The U.S. Government and The Netherlands Ministry of Defense are authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein. B.v.S. thanks Francesco Borsoi for fruitful discussions about latching.

CODE AVAILABILITY

The QArray+ Python package will be available on GitHub.

AUTHOR CONTRIBUTIONS

A.P. and B.v.S. developed the stochastic capacitance model. P.V. developed and implemented the spinless Hubbard model and the tunnel coupling in collaboration with B.v.S. E.D.N. and R.M. developed the charge truncation, supervised by P.V. and B.v.S. P.V. developed the open quantum system code under the supervision of T.P., B.K., N.A., and B.v.S. P.V. developed the multi-GPU acceleration under the supervision of T.P., B.K., and N.A. M.V. supervised A.P. and B.v.S. The manuscript was written by P.V. and B.v.S. with helpful contributions and comments from all authors.

15 DECLARATIONS

M.V. is a founding advisor of Groove Quantum BV and declares equity interests. N.A. declares an interest as a founder of QuantrolOx, which develops machine learningbased software for quantum control. All remaining authors declare no conflicts of interest.

[1] Hubert C. George, Mateusz T. Madzik, Eric M. Henry, Andrew J. Wagner, Mohammad M. Islam, Felix Borjans, Elliot J. Connors, J. Corrigan, Matthew Curry, Michael K. Harper, Daniel Keith, Lester Lampert, Florian Luthi, Fahd A. Mohiyaddin, Sandra Murcia, Rohit Nair, Rambert Nahm, Aditi Nethwewala, Samuel Neyens, Bishnu Patra, Roy D. Raharjo, Carly Rogan, Rostyslav Savytskyy, Thomas F. Watson, Josh Ziegler, Otto K. Zietz, Stefano Pellerano, Ravi Pillarisetty, Nathaniel C. Bishop, Stephanie A. Bojarski, Jeanette Roberts, and James S. Clarke. 12-spin-qubit arrays fabricated on a 300 mm semiconductor manufacturing line. Nano Letters, 25(2):793–799, December 2024. [2] Paul Steinacker, Nard Dumoulin Stuyck, Wee Han Lim, Tuomo Tanttu, MengKe Feng, Santiago Serrano, Andreas Nickl, Marco Candido, Jesus D Cifuentes, Ensar Vahapoglu, et al. Industry-compatible silicon spin-qubit unit cells exceeding 99% fidelity. Nature, 646(8083):81– 87, 2025. [3] Andreas Nickl, Nard Dumoulin Stuyck, Paul Steinacker, Jesus D. Cifuentes, Santiago Serrano, MengKe Feng, Ensar Vahapoglu, Fay E. Hudson, Kok Wai Chan, Stefan Kubicek, Julien Jussot, Yann Canvel, Sofie Beyne, Yosuke Shimura, Roger Loo, Clement Godfrin, Bart Raes, Sylvain Baudot, Danny Wan, Arne Laucht, ChihHwan Yang, Wee Han Lim, Andre Saraiva, Christopher C. Escott, Kristiaan De Greve, Andrew S. Dzurak, and Tuomo Tanttu. Eight-qubit operation of a 300 mm simos foundry-fabricated device, 2025. [4] Jacob F. Chittock-Wood, Ross C. C. Leon, Michael A. Fogarty, Tara Murphy, Sofia M. Patomäki, Giovanni A. Oakes, James Williams, Felix-Ekkehard von Horstig, Nathan Johnson, Julien Jussot, Stefan Kubicek, Bogdan Govoreanu, David F. Wise, M. Fernando Gonzalez-Zalba, and John J. L. Morton. Radio-frequency cascade readout of coupled spin qubits fabricated using a 300 mm wafer process, 2025. [5] Alexander Lidiak, Jacob Swain, David L. Craig, Joseph Hickie, Yikai Yang, Federico Fedele, Jaime Saez-Mollejo, Andrea Ballabio, Daniel Chrastina, Giovanni Isella, Georgios Katsaros, Dominic T. Lennon, Vincent P. Michal, Erik M. Gauger, and Natalia Ares. Virtual gates enabled by digital surrogate of quantum dot devices, 2025. [6] Rahul Marchand, Lucas Schorling, Cornelius Carlsson, Jonas Schuff, Barnaby van Straaten, Taylor L. Patti, Federico Fedele, Joshua Ziegler, Parth Girdhar, Pranav Vaidhyanathan, and Natalia Ares. End-to-end analysis of charge stability diagrams with transformers, 2025. [7] Anantha S. Rao, Donovan Buterakos, Barnaby van Straaten, Valentin John, Cécile X. Yu, Stefan D. Oost-

erhout, Lucas Stehouwer, Giordano Scappucci, Menno Veldhorst, Francesco Borsoi, and Justyna P. Zwolak. Modular autonomous virtualization system for twodimensional semiconductor quantum dot arrays. Physical Review X, 15(2), April 2025. [8] Cornelius Carlsson, Jaime Saez-Mollejo, Federico Fedele, Stefano Calcaterra, Daniel Chrastina, Giovanni Isella, Georgios Katsaros, and Natalia Ares. Automated all-rf tuning for spin qubit readout and control, 2025. [9] Justyna P. Zwolak, Thomas McJunkin, Sandesh S. Kalantre, J.P. Dodson, E.R. MacQuarrie, D.E. Savage, M.G. Lagally, S.N. Coppersmith, Mark A. Eriksson, and Jacob M. Taylor. Autotuning of double-dot devices in situ with machine learning. Physical Review Applied, 13(3), March 2020. [10] Donovan L Buterakos, Sandesh S Kalantre, Joshua Ziegler, Jacob M Taylor, and Justyna P Zwolak. Qdflow: A python package for physics simulations of quantum dot devices. SciPost Physics Codebases, page 065, 2026. [11] H. Moon, D. T. Lennon, J. Kirkpatrick, N. M. van Esbroeck, L. C. Camenzind, Liuqi Yu, F. Vigneau, D. M. Zumbühl, G. A. D. Briggs, M. A. Osborne, D. Sejdinovic, E. A. Laird, and N. Ares. Machine learning enables completely automatic tuning of a quantum device faster than human experts. Nature Communications, 11(1), August 2020. [12] Sandesh S. Kalantre, Justyna P. Zwolak, Stephen Ragole, Xingyao Wu, Neil M. Zimmerman, M. D. Stewart, and Jacob M. Taylor. Machine learning techniques for state recognition and auto-tuning in quantum dots. npj Quantum Information, 5(1), January 2019. [13] Tyler J. Kovach, Daniel Schug, M. A. Wolfe, E. R. MacQuarrie, Patrick J. Walsh, Owen M. Eskandari, Jared Benson, Mark Friesen, M. A. Eriksson, and Justyna P. Zwolak. Bootstrapping, autonomous testing, and initialization system for si/sige multi-quantum dot devices, 2025. [14] Justyna P. Zwolak, Jacob M. Taylor, Reed W. Andrews, Jared Benson, Garnett W. Bryant, Donovan Buterakos, Anasua Chatterjee, Sankar Das Sarma, Mark A. Eriksson, Eliška Greplová, Michael J. Gullans, Fabian Hader, Tyler J. Kovach, Pranav S. Mundada, Mick Ramsey, Torbjørn Rasmussen, Brandon Severin, Anthony Sigillito, Brennan Undseth, and Brian Weber. Data needs and challenges for quantum dot devices automation. npj Quantum Information, 10(1), October 2024. [15] Barnaby van Straaten, Joseph Hickie, Lucas Schorling, Jonas Schuff, Federico Fedele, and Natalia Ares. Qarray: A gpu-accelerated constant capacitance model simulator for large quantum dot arrays. SciPost Physics Codebases, October 2024.

16 [16] Tara Murphy, Katarina Brlec, Giovanni Oakes, Lorenzo Peri, Henning Sirringhaus, Henry Moss, M. Fernando Gonzalez Zalba, and David Wise. Rf-squad: A radiofrequency simulator for quantum dot arrays, 2025. [17] Jan Krzywda, Weikun Liu, Evert van Nieuwenburg, and Oswin Krause. Qdarts: A quantum dot array transition simulator for finding charge transitions in the presence of finite tunnel couplings, non-constant charging energies and sensor dots. SciPost Physics Codebases, January 2025. [18] Valentina Gualtieri, Charles Renshaw-Whitman, Vinicius Hernandes, and Eliska Greplova. Qdsim: A userfriendly toolbox for simulating large-scale quantum dot devices. SciPost Physics Codebases, January 2025. [19] Lucas Schorling, Pranav Vaidhyanathan, Jonas Schuff, Miguel J Carballido, Dominik Zumbühl, Gerard Milburn, Florian Marquardt, Jakob Foerster, Michael A Osborne, and Natalia Ares. Meta-learning characteristics and dynamics of quantum systems. arXiv preprint arXiv:2503.10492, 2025. [20] Jonas Schuff, Miguel J Carballido, Madeleine Kotzagiannidis, Juan Carlos Calvo, Marco Caselli, Jacob Rawling, David L Craig, Barnaby van Straaten, Brandon Severin, Federico Fedele, et al. Fully autonomous tuning of a spin qubit. Nature Electronics 9, pages 304–313, 2026. [21] Francesco Borsoi, Nico W Hendrickx, Valentin John, Marcel Meyer, Sayr Motz, Floor Van Riggelen, Amir Sammak, Sander L De Snoo, Giordano Scappucci, and Menno Veldhorst. Shared control of a 16 semiconductor quantum dot crossbar array. Nature Nanotechnology, 19(1):21–27, 2024. [22] Valentin John, Cécile X. Yu, Barnaby van Straaten, Esteban A. Rodrı́guez-Mena, Mauricio Rodrı́guez, Stefan D. Oosterhout, Lucas E. A. Stehouwer, Giordano Scappucci, Maximilian Rimbach-Russ, Stefano Bosco, Francesco Borsoi, Yann-Michel Niquet, and Menno Veldhorst. Robust and localised control of a 10-spin qubit array in germanium. Nature Communications, 16(1):10560, 2025. [23] Maxim De Smet, Yuta Matsumoto, Anne-Marije J. Zwerver, Larysa Tryputen, Sander L. de Snoo, Sergey V. Amitonov, Sam R. Katiraee-Far, Amir Sammak, Nodar Samkharadze, Önder Gül, Rick N. M. Wasserman, Eliška Greplová, Maximilian Rimbach-Russ, Giordano Scappucci, and Lieven M. K. Vandersypen. High-fidelity single-spin shuttling in silicon. Nature Nanotechnology, 20(7):866–872, June 2025. [24] C. H. Yang, A. Rossi, N. S. Lai, R. Leon, W. H. Lim, and A. S. Dzurak. Charge state hysteresis in semiconductor quantum dots. Applied Physics Letters, 105(18), November 2014. [25] Fabian Hader, Sarah Fleitmann, Jan Vogelbruch, Lotte Geck, and Stefan van Waasen. Simulation of charge stability diagrams for automated tuning solutions (simcats). IEEE Transactions on Quantum Engineering, 5:1–14, 2024. [26] W. G. van der Wiel, S. De Franceschi, J. M. Elzerman, T. Fujisawa, S. Tarucha, and L. P. Kouwenhoven. Electron transport through double quantum dots. Rev. Mod. Phys., 75:1–22, Dec 2002. [27] R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. K. Vandersypen. Spins in few-electron quantum dots. Rev. Mod. Phys., 79:1217–1265, Oct 2007.

[28] Floor van Riggelen-Doelman, Chien-An Wang, Sander L. de Snoo, William I. L. Lawrie, Nico W. Hendrickx, Maximilian Rimbach-Russ, Amir Sammak, Giordano Scappucci, Corentin Déprez, and Menno Veldhorst. Coherent spin qubit shuttling through germanium quantum dots. Nature Communications, 15(1), July 2024. [29] Andrew J Daley. Quantum trajectories and open manybody quantum systems. Advances in Physics, 63(2):77– 149, 2014. [30] Martin B Plenio and Peter L Knight. The quantum-jump approach to dissipative dynamics in quantum optics. Reviews of Modern Physics, 70(1):101, 1998. [31] Heinz-Peter Breuer and Francesco Petruccione. The theory of open quantum systems. OUP Oxford, 2002.

17 and (ii) using multi-GPU/multi-node execution for highthroughput dataset generation in regimes with strong interdot connectivity.

Appendix B: Comparison of simulators with QArray latching

FIG. 9. JIT execution time vs. tunnel-coupling density for a 20-dot array (64×64 grid, multi-GPU). The JAX just-in-time (JIT) execution time increases as the tunnelcoupling density grows, consistent with a larger effective connectivity that increases the number of tunnel-coupled transitions (and, consequently, the amount of work per voltage point) even when distributing the computation across multiple GPUs.

Appendix A: Computation complexity

Here we discuss the complexity of computing an N ×N charge stability diagram using the stochastic capacitance model for a system consisting of ndot dots. Per subinterval, there are 2ndot loading/unloading transitions and ndot (ndot −1) interdot transitions, i.e., ndot (ndot +1) transitions per sub-interval and hence Nr ndot (ndot + 1) candidate charge transitions per pixel. For each of these O(Nr n2dot ) transitions it is necessary to evaluate the electrostatic energy according to Equation 1, which requires multiplying an ndot × ndot matrix by a vector
of length 2 ndot , an operation with complexity O n . The perdot
pixel cost is therefore O Nr n4dot , and the overall computational complexity of the charge stability diagram is
O N 2 Nr n4dot . The computational complexity of this algorithm therefore scales more favorably with the number of dots than that of steady-state methods. We note that ⃗ caching C−1 dd Q at each pixel reduces the cost of evaluating the energy change of a single transition to O(ndot
), which would tighten the per-pixel scaling to O Nr n3dot . This polynomial scaling represents a dramatic improvement over steady-state solvers, which typically must evaluate the electrostatic energy of at least 2 ndot distinct charge states. Therefore, their complexity scales as at least O(2 ndot ). As shown in Fig. 9, JIT execution time rises rapidly at high tunnel-coupling density. In practice, dense tunnel graphs increase the number of nonzero couplings tij (and associated rates) that must be evaluated at each voltage point, which can dominate runtime even under multi-GPU parallelism. This benchmark highlights the importance of (i) exploiting JAX/XLA fusion

As per Table I, the original QArray [15] and QDFlow [10] have the ability to simulate latching in an unphysical fashion. They neglect the gate-voltage dependence of the relaxation rate and the possibility of relaxation through intermediate charge states. In Figure 10, we simulate the same double quantum dot (DQD) system as in Figure 1, though here the original QArray is used to simulate charge latching. The system is configured with tunneling rates ⃗Γτ = (1, 0), such that dot 1 readily exchanges charges with the reservoir while dot 2 remains isolated. Consider the latching that occurs when the (1, 1) → (1, 0) transition is hindered by a negligible loading/unloading rate. In the standard QArray simulation, the system remains trapped in the (1, 1) state until it crosses the (1, 0) → (0, 0) unloading boundary. In contrast, the more physically grounded QArray+ simulation remains in (1, 1) only until reaching the extended (1, 1) → (0, 1) unloading line. At this point, a rapid two-step transition (1, 1) → (0, 1) → (1, 0) becomes possible, facilitated by the fast loading/unloading of Dot 1 and the underlying inter-dot transition rates. In Figure 11 we simulate the same double dot system, but now with ⃗Γτ = (0, 0.1). With these rates, the QArray latching simulation breaks down entirely: in Figure 11(a), the charge state remains latched within narrow horizontal bands that streak across the entire sweep, bearing little resemblance to the extended latched regions produced by the physically grounded models in panels (b) and (d). In the QArray+ simulation, by contrast, the vanishing rate of dot 1 prevents the (2, 1) → (1, 1) transition from occurring, so instead there is an extended interdot transition corresponding to (2, 1) → (1, 2) → (1, 1).

Appendix C: Discrete-Time Steady-State Occupation

In this appendix, we derive the discrete-time steadystate occupation probability utilized in the latching model. We consider a single quantum dot coupled to one electronic reservoir and restrict our attention to two charge states, ⃗n0 and ⃗n1 . The dimensionless electrostatic energy change associated with a loading event is defined as x ≡ β∆U,

β = (kB T )−1 .

(C1)

18

42 ms

V1 (V)

QArray+ Lindblad

0.0

801 ms

-2.0 0.0 -2.0

V1 (V)

0.0

FIG. 10. Double-dot charge stability diagrams for reservoir couplings ⃗ Γτ = (1, 0), computed using four different methods. (a) The original QArray latching implementation, and (b)–(d) the QArray+ stochastic, spinless Hubbard, and Lindblad (open-system) models, respectively, for the same double-dot system. Dot 1 exchanges charge rapidly with the reservoir while dot 2 is isolated. Compute times for each 200 × 200 scan are indicated in the top-right corners.

1.

QArray

(b) 141 ms 0.0

Transition Rates and Probabilities

-2.0 0.0 -2.0 (d) QArray+ Hubbard 226 ms 0.0

V1 (V)

-2.0 -2.0

V1 (V)

(C2a) (C2b)

where f (x) = (1 + ex )−1 is the Fermi function and Γ denotes the bare tunnel rate. To model the system over discrete integration time steps of duration τ , we define the probabilities A and B, which represent the probability of the system remaining in charge state 0 and 1, respectively, during the interval τ: (C3a) (C3b)

V1 (V) QArray+ Lindblad

0.0

V1 (V)

0.0

199 ms

-2.0 0.0 -2.0

Steady-State Derivation

The discrete-time evolution over one integration step is described by the stochastic transition matrix T:   A 1−B T= . (C4) 1−A B The stationary occupation probability π1 satisfies the eigenvalue condition π = Tπ, where π = (1 − π1 , π1 )T . This leads to the balance equation: π1 = (1 − A)(1 − π1 ) + Bπ1 .

(C5)

Solving for π1 yields the steady-state occupation: π1 (x) =

A = exp [−Γτ f (x)] , B = exp [−Γτ f (−x)] .

44 ms

FIG. 11. Double-dot charge stability diagrams for reservoir couplings ⃗ Γτ = (0, 0.1), computed using four different methods. (a) The original QArray latching implementation, and (b)–(d) the QArray+ stochastic, spinless Hubbard, and Lindblad (open-system) models, respectively. Here dot 1 is fully decoupled from its reservoir and dot 2 is weakly coupled. Compute times for each 200 × 200 scan are indicated in the top-right corners.

2.

Electron exchange with the reservoir is governed by Fermi-Dirac statistics. The loading (Γ0→1 ) and unloading (Γ1→0 ) rates are given by Γ0→1 = Γf (x), Γ1→0 = Γf (−x),

-2.0 -2.0 (c) 0.0

QArray+ Stochastic

V2 (V)

V1 (V)

V2 (V)

-2.0 -2.0

(a) 0.0

V2 (V)

QArray+ Stochastic

V2 (V)

-2.0 0.0 -2.0 (d) QArray+ Hubbard 384 ms 0.0

V1 (V)

V2 (V)

-2.0 -2.0 (c) 0.0

(b) 144 ms 0.0

V2 (V)

V2 (V)

QArray

V2 (V)

(a) 0.0

1−A . 2−A−B

(C6)

The behavior of the system is fundamentally governed by the dimensionless “bandwidth” parameter Γτ , which determines the deviation from the continuous-time limit.

19

π1 (x)

0.8 0.6

35 30

300

25

310

20

320

0.4

15

0.2

10

0.0

sensor signal (a.u.)

Γτ = 10 0 Γτ = 10 1 Γτ = 10 2 Γτ = 10 3 Γτ = 10 4

vP2 (mV)

1.0

330 340 15

30

20

10

0

10

x = β∆U (dimensionless)

20

30

FIG. 12. Steady-state occupation π1 as a function of the dimensionless electrostatic energy x. The curves illustrate how the transition broadens and flattens with the parameter Γτ .

20

25

30 vP3 (mV)

35

40

FIG. 13. Experimentally measured CSD showing a double transition line, adapted from Ref. [28].

Appendix D: Experimentally measured CSD

Figure 13 shows an experimentally measured charge stability diagram exhibiting a double transition line at a dot–reservoir charge transition. This is the same feature reproduced by the stochastic capacitance model at finite temperature and fast lead tunneling in Figure 5(b), where resolving multiple stochastic transitions per pixel (Nr = 10) yields the time-averaged “double-step” signal discussed in section III.

Related documents

• QArray+: A physics-informed GPU-accelerated simulator for quantum dot arrays#924872
• QArray+: A physics-informed GPU-accelerated simulator for quantum dot arrays#920638
• Study High-Performance Computing Techniques for Optimizing and Accelerating AI Algorithms Using Quantum Computing and Specialized Hardware#84793
• Qiskit Machine Learning: an open-source library for quantum machine learning tasks at scale on quantum hardware and classical simulators.#489519
• Efficient Trotter-Suzuki Schemes for Long-time Quantum Dynamics#1128549