A hybrid MPS–HEOM method yields the minimum molecule count NT for disordered molecular polaritons to reach the thermodynamic limit, showing phonon timescales govern dark-state activation

The authors build an MPS–HEOM hybrid framework that represents the HEOM auxiliary density operators as a matrix product state, capturing non-Markovian vibrational relaxation together with Markovian cavity loss and external driving in one tensor-network framework, and simulate up to roughly 100 two-level systems in the intermediate-coupling regime with computational cost scaling linearly with system size. Previously, nonperturbative HEOM scaled exponentially and was restricted to the few-emitter regime, while methods reaching the thermodynamic limit often neglected dissipation and pumping; this work combines both, enabling numerically exact dynamics of disordered TC and HTC models from few emitters to the macroscopic limit within a single framework. The text reports a tensor-network structure with 2+N(K+1) sites, TDVP time evolution with GPU acceleration, about 60x speedup, roughly 0.3 s per timestep for the TC model and 2 s for the HTC model at N=10, plus convergence checks of HEOM truncation parameters (Nb=4, L=5 for dynamic disorder; Nb=8 for static disorder).

The authors introduce the convergence scale NT, the number of molecules required for photonic dynamics to reach the thermodynamic limit, determined by the time-normalized RMSE between consecutive TLS numbers with a threshold of 10^-4, and use it to quantitatively answer what minimum system size is needed for collective polaritonic systems to reach the thermodynamic limit. Experiments use macroscopic samples (generally >10^5 molecules) while most theoretical methods apply to systems with <20 molecules, leaving a gap without a quantitative yardstick; NT provides an operational criterion bridging the few-emitter and macroscopic regimes. Using the experimentally measurable average photon number ⟨a†a⟩ as the observable, with parameters aligned to BODIPY-Br-like organic polaritons (ωc=ω0=2.0, κ=Γ↓=0.02, Q=100), the work obtains NT=3 in the absence of disorder and systematically scans disorder strength σ and bath characteristic frequency γ.

The authors find that light–matter coupling disorder requires larger NT than frequency disorder, and that dynamic disorder is generally more demanding than static disorder; NT increases and then decreases as the bath becomes more Markovian (larger γ), a non-monotonic behavior reminiscent of a Kramers turnover, with the static (inhomogeneous) limit providing a lower bound for NT in this parameter range. Prior work lacked a quantitative account of how disorder affects convergence to the thermodynamic limit; this work links NT directly to disorder type, disorder strength, and phonon timescale, showing that phonon timescales control both the breakdown of collective behavior and the growth of NT. Based on the scans of NT versus σ and γ in Fig. 2 and the time dependence of NT in Fig. 3(d): at γ=0.08, NT first rises, peaks near t≈40, then decreases and saturates, whereas under static disorder and γ=0.02, NT is nearly time-independent.

The authors attribute these trends to disorder suppressing collective light–matter dynamics: frequency disorder couples bright and dark manifolds through the matter Hamiltonian, whereas coupling disorder directly breaks collective symmetry so the cavity mode couples to dark states, which become optically active gray states; the turnover in dark-state population is explained by second-order perturbative Fermi's golden rule, with the transfer rate set primarily by the spectral weight of the noise spectrum at the Rabi splitting. The work identifies the microscopic mechanism by which disorder activates non-collective degrees of freedom and establishes the suppression of collective behavior as the key mechanism governing thermodynamic convergence in disordered light–matter systems, highlighting the importance of polariton–vibrational timescale matching. The Supplemental Material derives the bright–dark coupling matrix element and the FGR rate, giving k_B→D=(N−1)/N·J(Ω_R)coth(βΩ_R/2), approaching J(Ω_R)coth(βΩ_R/2) in the thermodynamic limit; the perturbative theory reproduces the turnover seen in numerically exact results but underestimates the transfer rate in the small-γ regime.

The results are aimed at theorists and ab initio simulators studying collective molecular polariton behavior: within the TC/HTC models, a Debye-Drude vibrational spectral density, single-excitation initial conditions, and the chosen cavity loss/pumping rates, NT serves as a quantitative basis for choosing simulation system size and suggests that phonon engineering can selectively enhance dissipative pathways funneling excitations into the dark-state manifold. For experimentalists, it offers a framework for judging how disorder and phonon timescales affect when collective behavior breaks down.

This load is an incomplete reading: the specific numerical curves of Figs. 2 and 3 and Supplemental Figs. S1–S4 are known only through textual description, so figure details cannot be verified; the quantitative values of NT depend on the convergence threshold of 10^-4, the observation window tmax, and the chosen HEOM truncation parameters, so numbers may differ under other settings; the perturbative FGR underestimates the transfer rate in the small-γ regime, meaning quantitative interpretation there rests on the numerical results; and the conclusions are built on a specific spectral density form and parameters aligned to a single experimental system, so applicability to other molecular systems and multimode cavities remains an open question.