Abstract
Recent advances in numerically exact quantum dynamics methods have brought the dream of accurately modeling the dynamics of chemically complex open systems within reach. Path-integral-based methods, hierarchical equations of motion, and quantum analog simulators all require the spectral density (SD) of the environment to describe its effect on the system. Here, we focus on the decoherence dynamics of electronically excited species in solution in the common case where nonradiative electronic relaxation dominates and is much slower than electronic dephasing. We show that the computed relaxation rate is highly sensitive to the choice of SD representation─such as the Drude-Lorentz or Brownian modes─or strategy used to capture the main SD features, even when early-time dephasing dynamics remains robust. The key reason is that electronic relaxation is dominated by the resonant contribution from the high-frequency tails of the SD, which are orders of magnitude weaker than the main features of the SD and can vary significantly between strategies. This finding highlights an important, yet overlooked, numerical challenge: obtaining an accurate SD requires capturing its structure over several orders of magnitude to ensure correct decoherence dynamics at both early and late times. To address this, we provide a simple transformation that recovers the correct relaxation rates in quantum simulations constrained by algorithmic or physical limitations on the shape of the SD. Our findings enable a comparison of different numerically exact simulation methods and expand the capabilities of analog simulations of open quantum dynamics.