For molecular polaritons, disorder and phonon timescales control the activation of dark states in the thermodynamic limit

Collective light-matter systems host an extensive manifold of dark states whose role in the emergence of thermodynamic behavior remains poorly understood, especially in the presence of disorder and structured environments. Here, we develop a hybrid matrix product state-hierarchical equations of motion (MPS-HEOM) approach that enables numerically exact simulations of polariton dynamics from a few emitters to the thermodynamic limit under both static and dynamic disorder. This allows us, for the first time, to provide a quantitative and operational answer to the long-standing question of what is the minimum system size required to reach the thermodynamic limit in collective polaritonic systems. By introducing a convergence scale, N_T, i.e., the number of molecules required for the photonic dynamics to reach the thermodynamic limit, we show that dynamic disorder generally poses a greater computational challenge than static disorder. We attribute this behavior to the suppression of collective light-matter dynamics by disorder, which dynamically activates non-collective degrees of freedom. We further find that N_T exhibits a turnover behavior as the bath becomes more Markovian, as the bath timescales regulate bright-to-dark energy transfer and the involvement of dark and gray states. Hence, phonon timescales control both the breakdown of collective behavior and the growth of N_T. Our results establish the suppression of collective behavior as the key mechanism governing thermodynamic convergence in disordered light-matter systems.