Beyond the Quantum Regression Theorem in Variational Polaron Master Equations with Low-Dimensional Baths

While the quantum regression theorem (QRT) is the standard tool for computing multi-time correlation functions in open quantum systems, it relies on system-bath separability and an environment that remains in equilibrium, assumptions that are violated once dynamical correlations develop. Using the projection operator formalism, we derive an extension to the QRT that explicitly incorporates these correlation-induced corrections. We apply this framework to the variational polaron master equation for the spin-boson model in ohmic and super-ohmic regimes, where the polaron transformation mixes system-bath degrees of freedom to produce a non-thermal effective environment. Benchmarking against numerically exact tensor-network simulations demonstrates quantitative agreement for single- and two-time observables, including linear-response spectra, even at strong coupling. Our approach broadens the reach of analytic master equations to strong-coupling regimes, enabling treatment of multi-time observables where environmental memory effects and system-bath correlations are crucial.