Quantum and classical effects in system-bath correlations and optical line shapes

System-bath correlations play a pivotal role in the study of open quantum systems and have emerged as an important feature in distinguishing between quantum and classical behavior in energy transfer and spectroscopy. A variety of classical models have been constructed to test the limits of purely classical mechanics to capture such effects and, thus, by elimination to identify features that are unique to quantum systems. Yet numerical modeling in such classical frameworks has been largely limited to perturbative approximations, leaving some doubt regarding the extent to which these models are capable of capturing system-bath correlation effects. To overcome this gap, we here apply the numerically exact hierarchical equations of motion (HEOM) framework to examine dynamics under a ``classical density matrix'' framework that follows strictly Newtonian dynamics. Under the assumption of a Drude-Lorenz spectral density, we show that the resulting HEOM framework becomes identical to the stochastic Liouville equation. By comparing numerically exact dynamical simulations under this model with various ``flavors'' of the quantum HEOM formalism and with experimental absorption and fluorescence data, we examine the various roles that system-bath interactions play in thermalization and spectroscopic line-shape functions. Although the numerically exact classical simulations do reveal transient system-bath correlations, these are generally weaker than their quantum counterparts and vanish as the system approaches thermal equilibrium to produce a factorized equipartition state. For experimentally relevant parameters, we find that the classical treatment provides a fairly accurate treatment of absorption spectroscopy but that loss of the quantum detailed balance condition leads to a (sometimes dramatic) failure in describing fluorescence spectra.