🔬 mesoscale physics

Real-axis least-squares bath discretization for real-time transport in quantum-dot Josephson junctions

This paper introduces a real-axis least-squares bath discretization method within the driven Liouville–von Neumann framework for superconducting reservoirs, which significantly reduces computational costs while accurately resolving transient coherence, nonlinear transport features, and damping-induced decay rates in quantum-dot Josephson junctions without the artificial smearing typical of previous finite-broadening approaches.

Original authors: Ruixin Zhou, Bing Dong, Yiyan Wang

Published 2026-09-25

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Ruixin Zhou, Bing Dong, Yiyan Wang

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the microscopic world of superconducting materials, electricity flows without any resistance, a phenomenon that relies on electrons pairing up to move in perfect unison. When scientists build tiny devices using these materials, often sandwiching a small island of matter between two superconducting wires, they create a laboratory for studying how quantum particles behave under extreme conditions. These setups are crucial for developing future quantum computers and ultra-sensitive sensors. However, simulating how these systems behave in real time is notoriously difficult. The superconducting wires act as vast reservoirs of energy, and to simulate them on a computer, researchers must approximate this infinite sea with a finite number of points. If they use too few points, the simulation eventually loops back on itself, creating false echoes that ruin the results. If they try to smooth out these loops by adding artificial friction, they blur the delicate quantum features they are trying to study, effectively erasing the very physics they want to understand.

A team of researchers at Shanghai Jiao Tong University has developed a new way to navigate this dilemma, allowing for much faster and more accurate simulations of these quantum devices. Their work focuses on a specific type of junction where a single quantum dot connects to superconducting leads. To understand the behavior of such a system, the researchers had to solve a problem of representation: how to replace the continuous, infinite energy spectrum of a superconductor with a manageable set of discrete points without losing the unique properties of the superconducting gap, the energy range where no particles can exist. Previous methods struggled to balance the need for high resolution at the edges of this gap with the need to suppress numerical errors over long periods of time.

The team introduced a method that treats the selection of these discrete points not as a fixed rule, but as an optimization problem. Instead of placing points evenly or using a standard mathematical curve, they used a computer algorithm to find the best possible arrangement of energy levels, their weights, and their individual damping rates. This approach, which they call a real-axis least-squares bath discretization, allows them to compress the simulation significantly. They found that by carefully tuning these parameters, they could achieve results that match high-resolution reference simulations while using far fewer computational resources. In their tests, they reduced the time required to simulate a specific transport curve from several hours down to just twenty minutes, a speedup of more than twenty times, without sacrificing the accuracy of the physical predictions.

The researchers validated their method by checking it against several established phenomena known to occur in these junctions. They successfully reproduced the static relationship between current and the phase difference across the junction, a fundamental property of superconductors. They also accurately captured the complex steps in current that appear when a voltage is applied, known as multiple Andreev reflections, which serve as a fingerprint of the junction's transmission spectrum. Furthermore, they simulated the behavior of the system after a sudden change in voltage, observing the oscillations of trapped particles that had been reported in previous studies. Their analysis confirmed that the decay of these oscillations follows a predictable scaling law, where the rate of decay is directly linked to the damping parameters chosen in the simulation. This allowed them to resolve decay rates that were previously too small to distinguish from the numerical noise of older methods.

One of the most significant findings was the ability to simulate the locking of the system's frequency to an external alternating current, a phenomenon known as Shapiro locking. This includes both integer and fractional steps, where the system locks into rhythms that are fractions of the driving frequency. The researchers demonstrated that their method could capture these locking behaviors even in regimes where other simulation techniques become computationally impossible due to the complexity of the numbers involved. They showed that the positions of these locking steps were exact and robust, while the amplitudes, though slightly affected by the numerical broadening, followed the expected physical trends. This capability opens the door to exploring complex, multi-terminal geometries and systems with irrational driving frequencies, which were previously out of reach for real-time simulations.

The study also clarified the limitations of their approach. While the method excels at capturing the positions of energy levels and the timing of oscillations, the absolute size of the currents in certain strong-drive scenarios remains a qualitative estimate rather than a precise number. This is because the method introduces a small amount of numerical broadening to keep the simulation stable, which slightly alters the peak heights of the current. However, the researchers showed that this effect is controlled and predictable, allowing them to extrapolate the results to what the true, un-broadened values would be. They also noted that for very long simulations, the density of the discrete points must increase as the damping decreases, a trade-off that cannot be avoided but is managed efficiently by their optimized setup.

By turning the discretization of the superconducting environment into a flexible, optimized process, this work provides a powerful new tool for the field. It allows scientists to simulate the transient, non-equilibrium behavior of superconducting nanostructures with a level of detail and speed that was previously unattainable. The method recovers established physical laws and reveals new dynamical behaviors, offering a clear path forward for designing and understanding the next generation of quantum devices. The researchers have made their data and code available, ensuring that the community can build upon this foundation to explore even more complex interactions in the quantum realm.

Technical Summary: Real-axis least-squares bath discretization for real-time transport in quantum-dot Josephson junctions

Problem Statement
Efficient real-time simulation of driven superconducting nanostructures is hindered by the representation of the superconducting leads. Three conflicting requirements create a bottleneck:

Recurrence: Finite discrete baths suffer from Poincaré recurrences, contaminating slow observables after a time $t_{rec} \approx \pi N_b/D$ . Suppressing this requires increasing the mode number $N_b$ linearly with simulation time $T$ , leading to a computational cost scaling near $T^4$ for dense correlation-matrix propagation.
Spectral Resolution: The BCS gap-edge singularity and subgap coherence (critical for Andreev physics) demand high spectral resolution precisely where standard discretization is most difficult.
Damping Trade-off: Global damping suppresses recurrences but broadens the gap edge and artificially shortens coherent lifetimes. A single global damping rate cannot simultaneously preserve long-time fidelity and steady-state accuracy. Existing methods (Floquet, NRG, NEGF with fixed Dynes broadening, HEOM) either fail for transients, are computationally expensive for long-time dynamics, or introduce a resolution floor that obscures decay rates below the broadening scale.

Methodology
The authors extend the Driven Liouville–von Neumann (DLvN) approach to superconducting reservoirs (Bogoliubov–de Gennes leads). The core innovation is treating bath discretization as a constrained optimization problem rather than a prescribed grid.

Formalism: The system is modeled as a noninteracting quantum dot coupled to BCS leads. The DLvN equation is adapted to the Nambu basis, ensuring the damping structure preserves particle-hole symmetry (assigning the same damping $\gamma_j$ to particle and hole components of each mode).
Gauge Discipline: To maintain consistency with a static reference state $C_{ref}$ , the time-dependent bias/drive is introduced via the hopping-phase (Peierls) gauge. This keeps bath orbitals frozen, avoiding the inconsistency of rotating terms on damped sites.
Real-Axis Least-Squares (LS) Framework: The hybridization self-energy of a continuous lead is approximated by a sum of $N_b$ damped BdG modes. The mode energies ( $\varepsilon_j$ ), spectral weights ( $w_j$ ), and mode-resolved dampings ( $\gamma_j$ ) are determined jointly by minimizing a weighted least-squares error against the target self-energy.
Four Design Constraints:
1. Gap-clean weighting: High penalty on spectral weight inside the gap to prevent parasitic ohmic shunts and protect Andreev bound states (ABS).
2. Per-mode $\gamma$ : In-band modes receive large damping (absorption), while gap-edge modes receive small damping (protection of long-lived physics).
3. Density floor: A dense, equally spaced core near the gap edge prevents "edge beat" artifacts caused by discrete modes outliving the main signal.
4. Normal-lead constraint: For normal leads, a dense core is retained in the transport window to pin Fermi occupations correctly, with compression applied only to the wings.

Key Contributions

Explicit Hamiltonian Framework: The method constructs an explicit Hamiltonian representation of the bath with optimized parameters, generalizing naturally to multiterminal normal-superconducting geometries and offering a route toward interacting impurity solvers.
Independent Calibration of Broadening: The damping parameter $\gamma$ is treated not as a physical parameter but as an independently calibrated numerical broadening. Its leading long-time effect on trapped-quasiparticle decay follows a parameter-free scaling law $\kappa \simeq w_{bath}\gamma_{eff}$ .
Joint Refinement Path: The framework defines a controllable refinement path $(\gamma \to 0, N_b \to \infty)$ , allowing decay rates below the resolution scale of previous finite-broadening calculations to be numerically resolved.
Computational Efficiency: By optimizing the bath, the method achieves a reduction in computational cost by more than an order of magnitude (e.g., reducing a 7.7-hour simulation to 20 minutes) while maintaining accuracy relative to high-resolution references.

Results and Validations
The framework was validated against established results and high-resolution references across four scenarios:

Static Current-Phase Relations (CPR): The method recovers continuum free-energy CPR benchmarks. While static CPR convergence is slower for real-axis LS than imaginary-axis variants, it is sufficient for dynamics.
Multiple Andreev Reflections (MAR): The method reproduces the subgap staircase, parity effects (odd-order enhancement), and dispersion relations. It matches a high-resolution Gaussian reference (110 modes) with only 32 modes (3–4% relative deviation).
Quasiparticle-Trapping Dynamics: The method reproduces the oscillation frequency of trapped quasiparticles after a bias quench. Crucially, it validates the scaling law $\kappa \simeq w_{bath}\gamma_{eff}$ and resolves decay rates ( $\kappa \approx 5.7 \times 10^{-4}$ ) below the resolution scale of previous finite-broadening calculations.
Shapiro Locking: The method captures integer and fractional Shapiro steps. It confirms selection rules (e.g., $\pi$ -periodicity for $m=2$ locks) and demonstrates capability in regimes where frequency-domain Floquet methods are computationally prohibitive (irrational bias points, high-denominator rationals).

Significance and Boundaries
The paper claims that this framework provides a robust, scalable route for real-time transport in superconducting junctions, overcoming the trade-off between recurrence suppression and lifetime fidelity. It enables the study of transient coherence and nonlinear transport across broad time and parameter ranges.

Stated Limitations:

Interactions: The current work is strictly limited to noninteracting dots ( $U=0$ ).
Finite-Gap Strong-Drive Amplitudes: While positions and selection rules for Shapiro steps are robust, the absolute amplitudes of strong-drive locking in the finite-gap regime deviate at the $O(1)$ level from frequency-domain Floquet steady states. These amplitudes are reported as qualitative.
Resolution Requirements: In the extreme small- $\kappa$ regime, the required mode density scales as $1/\gamma$ . Compression optimizes coefficients but cannot circumvent this fundamental resolution requirement.
High-Frequency Observables: The specific "density-floor" bath variant used for long-time spectroscopy is not suitable for high-frequency spectral moments due to far-tail anchors.

The authors conclude that the method successfully recovers established physics while enabling new capabilities in parameter mapping and transient dynamics, with the boundaries of its quantitative accuracy explicitly mapped for future development.

More like this