Optical spectra in the condensed phase: Capturing anharmonic and vibronic features using dynamic and static approaches
Simulating optical spectra in the condensed phase remains a challenge for theory due to the need to capture spectral signatures arising from anharmonicity and dynamical effects, such as vibronic progressions and their induced asymmetry. As such, numerous simulation methods have been developed that invoke different approximations and vary in their ability to capture different physical regimes. Here we use several models of chromophores in the condensed phase and ab initio molecular dynamics simulations to rigorously assess the applicability of methods to simulate optical absorption spectra. Specifically, we focus on the ensemble scheme, which can address anharmonic potential energy surfaces but relies on the applicability of extreme nuclear-electronic timescale separation; the Franck-Condon method, which includes dynamical effects but only at the harmonic level; as well as the recently introduced ensemble zero-temperature Franck-Condon approach, which straddles these limits. We also devote particular attention to the performance of methods derived from a cumulant expansion of the energy gap fluctuations and test the ability to approximate the requisite time correlation functions using classical dynamics with quantum correction factors. These results provide insights as to when these methods are applicable and able to capture the features of condensed phase spectra qualitatively and, in some cases, quantitatively across a range of regimes.
I. INTRODUCTION<!>II. THEORY<!>A. Static approaches<!>Vertical excitation energy<!>Ensemble approach<!>B. Dynamical Approaches<!>Ensemble Zero-Temperature Franck-Condon approach<!>The cumulant expansion<!>III. MODEL SYSTEMS FOR THE CHROMOPHORE<!>A. Brownian Oscillator Model<!>B. Generalized Brownian Oscillator Model<!>C. Morse oscillator<!>D. Phenolate<!>IV. MODEL FOR THE CONDENSED-PHASE ENVIRONMENT<!>V. RESULTS AND DISCUSSION<!>A. Generalized Brownian oscillator model<!>B. Anharmonic Morse Oscillator Model<!>C. Phenolate<!>VI. CONCLUSIONS
<p>Optical spectroscopy is one of the workhorses for identifying and characterizing molecules and materials, as well as their energy relaxation processes. An accurate description of linear absorption spectroscopy [1][2][3][4][5] of a material constitutes an essential first step in understanding its nonlinear spectroscopy. [6][7][8][9][10][11][12] Yet, despite its central importance, linear absorption spectroscopy in the condensed phase remains a challenge for theory because accurate methods are required both for describing the electronic surfaces and obtaining the spectra. Accurate modeling of electronic surfaces requires capturing interactions between the chromophore and its environment on its ground and excited electronic states. These interactions range from weak dispersion effects in nonpolar solvents to the stronger hydrogen bonding and dipole interactions present for charged chromophores in water. In particular, for many chromophores multiple solvation shells must be included in the electronic structure calculation to converge the excitation energies, [13][14][15][16] requiring large electronic structure calculations. Here we focus on the latter challenge of obtaining the absorption spectrum from the electronic surfaces, which is infeasible to do exactly for anharmonic condensed phase systems because it requires the vibronic eigenstates that describe a) These authors contributed equally b) Electronic mail: tmarkland@stanford.edu c) Electronic mail: cisborn@ucmerced.edu the coupling of nuclear dynamics with electronic transitions. This infeasibility has led to the development of a variety of approximate methods that utilize the vertical energy gap or its thermal distribution, estimate the vibronic eigenstates (usually within a harmonic approximation), or employ the full dynamics of the system.</p><p>For the simple case of a chromophore that interacts weakly with its environment, the peak position of the optical absorption spectrum is often obtained using a single electronic structure calculation of the vertical energy gap, corresponding to the ground to excited state transition, for a static structure, sometimes solvated using a continuum model. 4,[17][18][19][20][21][22] To mimic the effects of spectral broadening due to interactions with the environment and transitions from thermally accessible ground electronic state vibronic eigenstates to the multiplicity of excited electronic state vibronic eigenstates, the delta function spectrum corresponding to the computed vertical energy gap may then be dressed using a Gaussian or Lorentzian distribution, or a convolution of these (Voigt profile). 1 The choice of these distributions is motivated by the inhomogeneous (Gaussian) and homogeneous (Lorentzian) limits corresponding to slowly or rapidly evolving nuclear coordinates, respectively, 23 coupled weakly to a classical vertical energy gap obeying Gaussian statistics. 23 The ensemble method 13,16,[24][25][26][27] builds on this idea by instead approximating the spectrum using the distribution of vertical energy gaps generated from a thermal sampling of nuclear configurations on the ground state potential energy surface (PES). The ensemble method thus accounts for anharmonicity in the electronic surfaces and can incorporate the effects of specific solvent interactions on the excitation energies, albeit at the expense of assuming a strict separation of electronic and nuclear timescales. 6 It is important to note that, although molecular dynamics is frequently used to generate the nuclear configurations in the ensemble approach, the dynamics is only used as a sampling approach, i.e., the time-ordering of the configurations is not utilized, and hence one could instead sample the nuclear configurations using Monte Carlo. As such, the ensemble approach employs purely static input and does not capture the influence of nuclear dynamics on electronic relaxation that leads to features such as vibronic fine structure.</p><p>The Franck-Condon (FC) method [28][29][30][31][32][33][34][35][36][37][38][39][40] and related schemes 3,5,[41][42][43][44] also only require the static input of the ground and excited state frequencies of the optimized structures. These methods incorporate dynamical information by exploiting the vibronic eigenstates to formally shift the time dependence to frequency space, thereby capturing vibronic progressions in the spectrum. However, these methods use a harmonic approximation to the ground and excited state PESs to estimate the eigenstates and therefore poorly describe the spectrum for systems with strong anharmonicity. Because optimized structures of the ground and excited states are required for calculating the frequencies, explicit treatment of the solvation environment of the chromophore presents a challenge. 45 In particular, finding a minimum for an explicitly solvated chromophore can be difficult and in many cases impossible. Even if such a minimum can be found, the high dimensionality of the explicitly solvated system makes construction of the Hessian computationally costly and the significant geometry changes that occur in low frequency modes upon going from the electronic ground to the excited state often invalidates the harmonic approximation on which the FC expressions rely. This complication limits most FC calculations to a continuum description of the solvent where the minimum energy structures are well defined. Recently, the ensemble-zero-temperature Franck-Condon (E-ZTFC) method has been introduced which combines the appealing aspects of the ensemble and FC approaches. 22,46,47 Although originally introduced employing physically motivated arguments, E-ZTFC can be formally derived from the dynamical formulation of the ensemble and FC approaches, as is shown in Sec. II B 2, thereby elucidating the specific approximations that underlie this scheme.</p><p>The aforementioned methods do not require dynamical input. However, this information can be exploited using methods such as those based on a cumulant expansion with respect to electronic vertical energy gap fluctuations. 6,48 Methods derived from the cumulant expansion allow one to construct the absorption spectrum from the autocorrelation functions of the fluctuations of the energy gap operator. Because these methods employ time autocorrelation functions sampled on the full ground state PES, they account for both anharmonic effects and the dynamical effects that lead to vibronic spec-tral features. For many systems, a practical balance between expense and accuracy is to truncate the cumulant expansion at second order. The resulting second-order cumulant scheme is exact when the statistics of energy gap fluctuations are Gaussian and requires only the calculation of the equilibrium energy gap fluctuation (single time) autocorrelation function. 6,49 It should be noted that reliance of the second-order cumulant approach on the Gaussian nature of the fluctuations of the energy gap does not imply that the PESs need to be harmonic (see Sec. II B 3 and Refs. 48 and 50). The second-order cumulant scheme has been heavily exploited to reproduce and elucidate the vibrational spectroscopy of probes in the condensed phase [51][52][53][54] and has recently also found use in directly simulating electronic spectra in the condensed phase. [55][56][57] Beyond linear spectroscopy, it can also be used to investigate nonlinear spectroscopy and energy transfer in molecular systems. 6,49 However, the second-order cumulant has been observed to give poor results when compared to experimental vibrational spectra for several systems, including N-methylacetamide in methanol, 58 amide and carbonyl in water, 59 and the OH stretch in deuterated water. 60 With the advent of the use of the cumulant expansion for simulating linear optical absorption spectroscopy of condensed phase systems, it is timely to pursue a rigorous assessment of when the second-order cumulant approach can be expected to perform well, when including the third-order cumulant correction can improve the results, and how accurately the correlation functions required by these methods can be captured using classical dynamics with quantum correction factors.</p><p>Here we provide a rigorous assessment of the applicability of and links between a series of spectral simulation methods, including those utilizing dynamical input and those that require only static input, for the calculation of absorption spectra for systems in the condensed phase. By showing a unified derivation arising from a dynamical formulation, we elucidate the approximations underlying these approaches and their region of applicability, highlighting the physical systems for which these methods can be expected to perform well. To complement this analysis, we benchmark the performance of these methods on model systems chosen to provide stringent tests of the versatility of these methods, showcasing the capabilities of the second-and third-order cumulant expansion. In summary, this work provides a map for the advantages and disadvantages of a series of methods for modeling absorption spectra, specifies the approximations that underlie them, and demonstrates their ability to capture important physics that determines the shape of the absorption spectrum, including the anharmonicity of PESs and the dynamical effects that result in vibronic fine structure.</p><!><p>We are interested in predicting the linear absorption spectrum of a chromophore in solution where the spectrum is dominated by the transition from the ground to first excited electronic state. In this case, the Hamiltonian of the molecule and its environment takes the general form,</p><p>where g and e denote the electronic ground and first excited state, respectively, H g (q) and H e (q) are the nuclear Hamiltonians corresponding to the ground and excited electronic states, and V eg (q) = V ge (q) † are the transition dipole operators, which are generally dependent on the nuclear configuration.</p><p>For this system, the linear absorption spectrum can be written as,</p><p>where α(ω) is an ω-dependent prefactor whose form is dependent on the particular definition of the absorbance σ(ω) used. 6,43 For the purpose of this work, we use α(ω) = 1, such that the reported absorption spectra are directly proportional to the Fourier transform of χ(t).</p><p>The linear response function χ(t) is given by χ(t) = Tr ρ(0)e iHt V (q)e −iHt V (q) = Tr n ρ g (q, p)V ge (q, t)e</p><p>In the second line, we have performed the trace over electronic degrees of freedom, Tr e , leaving only the trace over nuclear degrees of freedom Tr n , moved into the interaction picture, and employed the Dyson expansion to express the response function as a nuclear average of the time-integrated energy gap operator</p><p>where V g (q) and V e (q) are the ground and excited state PESs and ω 0 eg is the excitation energy separating the minima on these two surfaces (also called the adiabatic energy gap). Upon moving to the interaction picture, the time dependence of an arbitrary (nuclear) operator O(p, q) is taken with respect to the nuclear Hamiltonian corresponding to the electronic ground state, O(t) = e i Ĥg( q,p)t O(p, q)e −i Ĥg( q,p)t .</p><p>For the absorption process, the initial condition with respect to which the average in Eq. ( 3) is performed is the direct product of the ground electronic state and the canonical density for the nuclear degrees of freedom on the ground PES,</p><p>The angled brackets in the last line of Eq. ( 3) denote the equilibrium average with respect to the nuclear degrees of freedom, i.e., ... ≡ Tr n [ρ g ...]. Equations ( 2) and ( 3) are the fundamental expressions required to calculate the optical absorption spectrum of a chromophore in its environment. In the following sections, we show how a number of commonly used approximate methods can be obtained from these expressions and elucidate the connections between them, emphasizing their ability to incorporate the effect of the environment on the excitation energy, treat anharmonic PESs, and include dynamical information, which permits the description of, for example, vibronic line progressions.</p><!><p>The approaches presented in this section rely on the assumption that the nuclear motions are slow in comparison to the timescale over which the absorption spectrum is taken. This limit is appropriate when the absorption spectrum is dominated by inhomogeneous broadening arising from static disorder, i.e., the main influence of the environment on the absorption spectrum is to provide the distribution of nuclear configurations and environments that modify the energy gap corresponding to the electronic transition. As such, the nuclear coordinates that appear in the response function in Eq. (3) are effectively frozen, i.e., q(t) ≈ q(0), which leads to the following simplification,</p><p>where</p><p>is the representation of the canonical nuclear operator in real space. In going from the first to second lines of Eq. ( 6), we have exploited the phase space representation of static averages. Substituting Eq. ( 6) into Eq. (2) then yields the following expression for the absorption spectrum,</p><p>(8) The second line of Eq. ( 6) represents a convenient starting point for a range of exact to approximate treatments.</p><p>For example, by using path integral molecular dynamics (PIMD) 61 simulations to sample the integral, one can obtain the exact quantum mechanical result for the response function in the inhomogeneous limit for a given level of electronic structure treatment of the ground and excited state PESs and vertical energy gap.</p><!><p>The simplest static approach can be obtained by setting f (R) = δ R − R 0 g in Eq. ( 7), which is exact in the classical limit at zero temperature, and where R 0 g corresponds to the energy minimum on the ground state PES. Upon substitution into Eq. ( 8), this yields</p><p>Thus, the absorption spectrum reduces to a δ-function placed at a vertical excitation energy computed for a single representative nuclear configuration R 0 g . To generate a spectrum with a finite width, δ(ω − U (R 0 g )) is often replaced by a Gaussian function, where the the standard deviation is chosen in a phenomenological way. 1 Because finding a single well-defined global minimum of the ground state PES is impractical in disordered condensed phase systems, R 0 g is often computed for a chromophore embedded in continuum solvent. 4,20,62 Hence, the vertical excitation energy approximation does not capture any spectral width beyond a phenomenological broadening parameter and typically does not account for any direct chromophore-environment interaction beyond a PCM. However, the approach requires only a single evaluation of the vertical energy gap, U (R 0 g ), making it the only tractable method if high-level wavefunction schemes are required to compute reliable excited state energies for the system of interest.</p><!><p>In the ensemble approach, one tries to obtain a finite temperature representation of the nuclear distribution in Eq. (8). For a given choice of the ground state PES, this can be obtained exactly using PIMD 61 or by sampling from the Wigner distribution. 63 However, due to the expense associated with generating the exact quantum nuclear distribution, classical MD is often used to sample the configurations, i.e., f (R) = dP e −H(P,R) Z which corresponds to the classical configurational distribution function.</p><p>Even when classical nuclear dynamics is employed, the need to perform sufficiently long ab initio molecular dynamics simulations on the ground state PES is challenging. Hence in practice these configurations are often generated using lower level methods such as force fields which leads to an inconsistency between the vertical energy gap evaluation and the underlying sampling. 55,64,65 The configurations used to generate the spectrum each require the calculation of the vertical energy gap and corresponding transition dipole. In practice, hundreds of excited state calculations are required to produce a converged spectrum. It is therefore common practice 46 to generate a continuous spectrum by broadening the δ-distribution in Eq. ( 8) to normal distributions centered around U (R) with a standard deviation chosen to provide a smooth spectrum.</p><!><p>The approaches described in this section, in contrast to those presented in Sec. II A, can be cast into a dynamical form where the time dependence of the nuclear coordinates is included, and hence go beyond the inhomogeneous limit.</p><p>1. Franck-Condon approach and its zero-temperature limit</p><p>The FC method [28][29][30][31][32][33][34][35][36][37][38][39][40] relies on the harmonic approximation to the ground and excited state PESs. In addition, this approach assumes the validity of the Condon approximation, i.e. that the transition dipole does not depend on the nuclear configuration, 66,67 V eg (q) ≈ V eg .</p><p>Several approaches exist that relax the Condon approximation by allowing V eg to vary linearly with nuclear degrees of freedom, 5,42-44 but these correction terms are mainly relevant for weak optical transitions and a detailed discussion is beyond the scope of this work. The response function for the FC method is</p><p>Here {|g ⊗ |ν g , |e ⊗ |ν e } are the vibronic eigenbases of the Hamiltonian, Eq. (1a), within the harmonic approximation such that H g (q) |ν g = E νg |ν g and H e (q) |ν e = E νe |ν e , Ω νe,νg = E νe − E νg is the energy gap between the ν g th and ν e th vibronic states on the ground and excited electronic states, respectively. We emphasize that this energy gap between vibronic states is not the same as the vertical excitation energy at nuclear configuration q, i.e., Ω νe,νg = U (q). Substituting Eq. ( 11) into Eq. (2) yields the FC spectrum</p><p>As in the case of the vertical excitation gap and ensemble methods, the delta functions in Eq. ( 12) are often broadened using Gaussian distributions.</p><p>For chemical systems where the frequencies characterizing its normal modes are larger than the thermal energy, one can instead calculate the spectrum in the zero-temperature limit. In this limit, the sum over initial states collapses to a single contribution from the ground vibronic state of the ground PES, i.e., ρ νg → δ 0g,νe , which leads to the following expression for the absorption spectrum (13)   and its associated zero-temperature FC (ZTFC) response function,</p><p>The FC approach has been shown to yield accurate absorption spectra in a variety of systems, 41,68-70 particularly in the case of small, rigid molecules in weakly interacting solvent environments, but can yield inaccurate spectra if the chromophore contains a number of anharmonic low-frequency modes that couple significantly to the electronic excitation. 71</p><!><p>The shortcomings of the FC approach, especially in condensed phase systems with strong solute-environment interactions, have motivated the recent development of a combined ensemble plus zero temperature Franck-Condon (E-ZTFC) approach. 22 The approach is aimed at retaining the main advantages of the ensemble approach, such as the explicit sampling of chromophoreenvironment interactions and low-frequency anharmonic motion of the chromophore coupled to its environment, while accounting for dynamic effects in the spectrum due to fast chromophore vibrations. In the following, we demonstrate how this method can be obtained starting from the dynamic formulation of the absorption spectrum.</p><p>To arrive at the E-ZTFC expression for the absorption spectrum, we begin from the expression for the response function, Eq. ( 3), and employ the ansatz that the response function, and therefore the spectrum, can be thought of as being composed of two largely noninteracting contributions. The first accounts for the inhomogeneous temperature-dependent limit and is embodied by the response function corresponding to the en-semble method, Eq. ( 8). The second component, consisting of the ZTFC response, Eq. ( 14), is temperatureindependent and accounts for the vibronic contribution to the spectrum. Hence, the E-ZTFC response function can be constructed as follows,</p><p>where the time-dependent function γ(t) is constrained such that the E-ZTFC response function, χ E−ZTFC (t), recovers either the ensemble or ZTFC response functions in the appropriate limits. Specifically, in the inhomogeneous limit where the nuclear degrees of freedom are approximately static over the timescale of the absorption process, χ E−ZTFC (t) must reduce to χ E (t), whereas in the zero-temperature harmonic limit, χ E−ZTFC (t) must reduce to χ ZTFC (t). Hence,</p><p>is the inhomogeneous limit of the ZTFC response function, Eq. ( 14), the form of which can also be obtained by taking the zero-temperature limit of the ensemble response function, Eq. ( 8). In the classical limit, i.e., when classical MD is used to sample the ground state PES in the ensemble method, Eq. ( 16) takes a particularly simple form,</p><p>where R 0 g is the nuclear configuration at the minimum of the ground state PES and U (R 0 g ) is the vertical excitation energy at this position.</p><p>To obtain the quantum limit, i.e., when using PIMD to sample the ground state PES in the ensemble method, Eq. ( 16), we exploit the fact that within the FC spectrum, the most general form of the parametrized Hamiltonian corresponds to the generalized Brownian oscillator model (GBOM), introduced in Sec. III B, for which we have provided response functions in Sec. II of the supporting information (SI). Thus, the quantum limit of Eq. ( 16) corresponds to the β → ∞ and t → 0 limits of the GBOM response function, Eq. ( 17) in SI Sec. II B, for the parametrized version of the GBOM arising from the normal mode analysis used in the FC approach.</p><p>Because the classical version of the ensemble approach is more widely used, we continue our discussion assuming the applicability of Eq. ( 17), but note that one can similarly obtain the quantum version by using the expression for χ inh ZTFC (t) provided in SI Sec. II C. Substituting Eqs. ( 15) and ( 16) into Eq. ( 2), we obtain,</p><p>which recapitulates the E-ZTFC expression for the absorption spectrum in Eq. ( 9) in Ref. 22. This derivation elucidates that the E-ZTFC method assumes that there are two major sources of broadening and structure in the absorption spectrum. The first corresponds to contribution arising from the ensemble method, which is assumed to capture all temperature related effects and is by construction static. The second corresponds to the ZTFC component, contains all relevant dynamical effects, and is assumed to contribute only at zero temperature.</p><!><p>The cumulant scheme 6,48 is a dynamics-based approach that re-expresses the response function within the Condon approximation as an infinite expansion with respect to cumulants (combinations of moments) of the energy gap fluctuation, δU (q, t) = U (q, t) − ω av eg ,</p><p>where</p><p>is the mth order approximation to the full lineshape function G ∞ (t). Here, g n (t) is the nth order cumulant and ω av eg = U (q) is the thermally averaged energy gap operator. To obtain explicit expressions for each cumulant, one expands both the time-ordered exponential on the first line of Eq. ( 19) and the exponentiated sum of cumulants in the second line and matches terms order by order in δU . When the energy gap fluctuation operator, δU (q), obeys Gaussian statistics, g n (t) = 0 for n ≥ 3. 6,49 The mth order cumulant approximation to the response function can then be expressed in terms of the mth order approximation to the lineshape function,</p><p>For the second-order cumulant, i.e., g 2 (t) = G 2 (t), where the expansion is truncated at second-order, the lineshape function takes the form,</p><p>The 2nd order lineshape function depends on the spectral density, J (ω), which can be expressed as follows</p><p>where θ(ω) is the Heaviside function and</p><p>is the energy gap fluctuation equilibrium time correlation function. However, the exact time correlation function C δU (t) can be calculated only for simple model systems. Although many schemes have been developed to approximate equilibrium time correlation functions, the standard approximation in the application of secondorder cumulant approach to general anharmonic systems is to use the classical version of C δU (t) to construct an approximate quantum analogue. 6,[72][73][74][75][76][77] Specifically, noting that the quantum Kubo-transformed correlation function obeys the same symmetries as classical equilibrium time correlation functions, i.e., both are real, even functions of time, the classical equilibrium time correlation function is used as an approximation to the quantum Kubotransformed analogue. 78,79 One can then employ general properties of quantum equilibrium time correlation functions to express the imaginary component of the energy gap fluctuation autocorrelation function, Im[C δU (t)], in terms of its Kubo transformed form, C kubo δU (t), leading to the following expression for the spectral density,</p><p>The specific choice of quantum correction factor employed in Eq. ( 25), while not unique, 73 has been widely used in the context of computing linear spectra of chromophores in complex environments, [55][56][57] and has shown good performance in constructing spectral densities for studying the energy transfer in multi-chromophore pigment-protein systems. [80][81][82][83][84] Hence, to calculate the linear absorption spectrum via the second-order cumulant approximation using the classical approximation for the energy gap fluctuation autocorrelation function, one calculates C cl δU (t) by calculating the vertical energy gap, U (q), over the course of a molecular dynamics simulation on the ground state PES and uses Eq. ( 25) to construct the lineshape function, Eq. ( 22), which can be used to construct the response function, Eq. ( 19). Finally, one Fourier transforms the computed response function, Eq. ( 19), to obtain the absorption spectrum, Eq. ( 2).</p><p>The second-order cumulant approximation therefore relies on two central approximations. First is the secondorder truncation of the cumulant expansion, which is exact for simple model systems like the Brownian oscillator model (BOM) introduced in Sec. III A. However, like Marcus theory, this approximation does not require that the underlying PESs be harmonic. Instead, the secondorder truncation implies that the time-dependent statistics of the energy gap fluctuation operator be Gaussian. One can, at least in principle, relax this approximation by incorporating higher-order cumulants. The second approximation consists of approximately constructing the quantum mechanical energy gap fluctuation autocorrelation function from its classical counterpart.</p><p>In the following, we calculate absorption spectra using both the second-and third-order cumulant approximations. The third-order truncation adds a correction to the response function,</p><p>, which is a functional of the two-time quantum correlation function of the energy gap operator δU , where the exact form can be found in SI Sec. I.</p><p>A major difficulty that the third-order cumulant presents is the need to calculate a multi-time quantum correlation function,</p><p>As in the case of the one-time correlation function, Eq. ( 24), needed to construct the second-order cumulant, calculation of Eq. ( 26) raises the question of how to best approximate multi-time quantum correlation functions in terms of their classical counterparts. Recently, an approximate link has been derived 85 to connect the Fourier transform of the double Kubo-transformed quantum correlation function to its classical counterpart,</p><p>where ω = ω+ω . Here we assess this recent development in the context of linear optical spectroscopy.</p><p>Although the uncontrolled approximations made in the cumulant approach in terms of constructing quantum correlation functions from their classical counterparts require careful considerations to test their applicability in realistic systems, the method is appealing because it can explicitly account for the coupling between the system and its condensed phase environment, in contrast to the FC approach. It also includes dynamic effects and does not suffer from the errors introduced the E-ZTFC approach due to decomposing the spectrum into temperature-independent dynamic and finitetemperature inhomogeneous effects. Thus, although the cumulant approach is the computationally costliest approach considered in this work (requiring the explicit propagation of the system on the ground state PES and the evaluation of vertical excitation energies at each timestep to evaluate the correlation functions required), it has the potential of providing improved spectra, especially in complex environments that are not well-represented by a classical continuum model.</p><!><p>In the following, we introduce several model systems for the chromophore degrees of freedom that are sufficiently simple as to allow for the exact (or very accurate) calculation their absorption spectra. The model systems include simple Hamiltonians where the ground and excited PESs can be decomposed into multidimensional harmonic surfaces, the Morse oscillator, and the phenolate ion (see Schematic in Fig. 1). Although we do not present spectroscopic results for the simplest model of a chromophore, the Brownian oscillator model (BOM), we describe this simple model so that it is clear how the generalized BOM differs from this model. Additionally, the condensed phase environment we describe in the next section is for a BOM-type environment.</p><!><p>The simplest of these models, the Brownian oscillator model (BOM) corresponds to nuclear Hamiltonians on the electronic ground and excited states consisting of noninteracting harmonic oscillators with shifted minima,</p><p>As defined at the beginning of Sec. II, ω 0 eg is the energy difference between the minima of the multidimensional ground and excited state PESs or adiabatic energy gap, pj and qj are the momentum and coordinate operators, ω j is the frequency, and K j is the shift of the minimum in going from the ground to the excited state PES of the jth oscillator. For this model system, both the FC approach within the harmonic approximation and the second-order truncation of the cumulant expansion are exact, i.e., Eq. ( 22), is the exact form of the lineshape function. Here, the energy gap fluctuation operator is linear in the oscillator coordinates, δU (q) = U (q) − ω av eg = j ω 2 j K j qj , where</p><p>Harmonic oscillators are sufficiently simple to allow analytical solutions to their quantum dynamics; the exact expression for the spectral density, J (ω), in Eq. ( 23) 6 is</p><p>Because the second-order cumulant expansion is exact for the nuclear Hamiltonian in Eqs. (28a)-(28b), it is common practice to interpret the spectral density as reporting on the frequency weighted difference of the shift in the minima of the Brownian oscillator model modes. Often the peaks in the spectral density are matched to normal modes of the chromophore. These comparisons suggest that the second-order cumulant is most applicable to harmonic systems of the Brownian oscillator model form and belie its generality and broad applicability to generally anharmonic systems. For this reason, in Sec. V, we test the applicability of the second-order cumulant approach to systems that go beyond the BOM.</p><!><p>The generalized form of the BOM (GBOM) also possesses harmonic ground and excited state PESs, but the surfaces may have different curvatures and the noninteracting harmonic modes on the electronic excited state  can be rotated with respect to the ground state,</p><p>where ω g,j and ω e,j are the frequencies of the ground and excited state normal modes and the coordinates of the normal modes on the ground and excited state PESs can be connected via a linear transformation called the Duschinsky rotation matrix, J, 86</p><p>The GBOM is the most general form of the Hamiltonian that results from the harmonic expansion of the ground and excited state PESs, and is therefore the output of the standard normal mode analysis of these surfaces performed in many electronic structure packages. This normal mode data provides the basis of computing linear spectra for complex molecules in the FC approach, 43,44,87,88 since it is possible to construct exact closed expressions for its linear response function using the harmonic oscillator propagator x|e iHt |y . The relatively simple form of the GBOM allows for a derivation of analytic expressions for all approaches to linear spectroscopy discussed in Sec. II (See SI Sec. II). In particular, the harmonic form of the PES means that exact expressions for the one-and two-time quantum correlation functions can be derived (SI Sec. II. D), allowing for a rigorous assessment of the errors introduced when constructing the cumulant expansion from classical correlation functions, as is commonly done for applications to realistic systems. 56,57 Although the spectrum of the GBOM can be solved analytically for an arbitrary number of modes, we limit the discussion here to a simple two-mode system. The two-mode system is the simplest system that still retains all the non-linearities in the energy gap operator introduced by both Duschinsky rotations and differences between ground and excited state normal mode frequencies. The ground state frequencies and displacement vector K of the GBOM used in this work are given in Table I. The excited state frequencies differ from the ground state frequencies by a factor ∆ω, with ω e,1 = ω g,1 (1 + ∆ω) and ω e,2 = ω g,2 (1 − ∆ω). The value of ∆ω is chosen to be 5% to match realistic normal mode frequency changes between the ground and excited state observed in simple chromophores. For the two-mode system, the Duschinsky rotation linking the ground and excited state normal mode coordinates is simply given by the 2D rotation matrix</p><p>and a value of φ = 10°is used for the GBOM results presented in this work. This choice is again informed by typical parameters for realistic systems, where some mode mixing occurs but the Duschinksy matrix is found to be diagonally dominant. 43</p><!><p>Thus far, we have focused exclusively on harmonic models. However, the PESs of real chemical systems are generally anharmonic. The anharmonic Morse oscillator potential takes the form,</p><p>where D is the well depth, α is the anharmonicity parameter (with the limit of α → 0 recovering the harmonic oscillator), and r eq is the equilibrium distance between the bonded atoms that minimizes the potential. For simplicity, in this work we consider only a single anharmonic mode. We fit the parameters for the Morse potentials describing the ground and excited state surface to a scan over the bond stretching mode of the carbon monoxide molecule. The calculations were performed in Gaussian 89 using density functional theory (DFT) to describe the ground state and time-dependent DFT (TDDFT) in the Tamm-Dancoff approximation 90 to describe S 1 excited state PES with the CAM-B3LYP 91 functional and the 6-31+G* 92 basis set. Further details of the fit and the resulting Morse potential parameters can be found in SI Sec. IV. We also consider a second, more anharmonic parameterization of the Morse oscillator, where α is increased by a factor of 1.8 from its value for CO for both ground and excited state PESs, and the parameters D are divided by a factor of 1.8 2 , guaranteeing that the Taylor expansion of the PESs around their respective minima remains unchanged up to second order. Due to the anharmonicity in the Morse potential, the FC approach based on a parameterized GBOM is no longer exact for this system. However, because the wavefunctions and energy levels of the Morse oscillator have analytic expressions, 93 it is still possible to obtain the exact response function of the system by evaluating Eq. ( 11) directly. In this work, we evaluate the linear response function for the Morse oscillator by computing the FC integrals of ground and excited state nuclear wavefunction overlaps numerically. A total number of 50 vibrational states on both energy surfaces are considered when computing the overlaps for the CO system. For the more anharmonic Morse oscillator test system, we instead include all bound vibrational states in the calculation of the FC overlaps.</p><p>We compare this exact spectrum to the other approaches. For the FC approach within the harmonic approximation, we use the expression for the Morse oscillator written as a harmonic oscillator with effective frequency</p><p>where µ is the effective mass of the vibrational mode.</p><p>To compute the spectra for this system in the cumulant approach, we performed classical MD on the onedimensional ground state PES using the i-PI code. 94,95 A trajectory of 400 ps, sampled every 2 fs, was used to construct the classical correlation functions. The quantum correlation functions were then constructed by applying the quantum correction factors introduced in Sec. II B 3. Details on the simulation parameters used in the i-PI code can be found in SI Sec. IV.</p><!><p>Lastly, we move away from model systems to study the S 0 → S 1 transition of the phenolate (C 6 H 5 O) − anion. Here, we compare the cumulant results obtained by performing dynamics on the ab-initio PES with the results obtained by parameterizing a GBOM, as is commonly done in FC approaches to the lineshape. To parameterize an effective GBOM for the FC spectrum within the harmonic approximation, we computed the ground and excited state optimized geometries and normal modes of phenolate in vacuum at the CAMB3LYP/6-31+G* level of theory using Gaussian, 89 where the Tamm-Dancoff approximation 90 is applied to the excited state TDDFT calculations. The shift vector K and the Duschinksy matrix J were also directly computed in Gaussian. To compute ensemble and cumulant spectra that consider the full PES of the system, without relying on a parameterization of the GBOM within the harmonic approximation, we performed an ab-initio MD simulation of phenolate in vacuum using the Terachem code. 96,97 The same basis set and DFT functional as for the normal mode calculations was used to perform the MD. Detailed simulation parameters for the ab-initio MD can be found in SI Sec. V. Three independent 50 ps trajectories were generated from the simulations. From these trajectories, the chromophore coordinates were extracted every 2 fs and the S 1 vertical excitation energy was calculated using the TeraChem code 98 in the Tamm-Dancoff approximation. The total of 75,000 data points was then used to construct either the ensemble spectrum or the one-and two-time classical correlation functions for the cumulant approach. Although the phenolate electronic structure and AIMD calculations are performed in vacuum, the simulated spectra include the effects of the BOM model of the environment described below.</p><!><p>Effects of the condensed phase environment in this work are introduced via immersion of the model chromophore system into a BOM-type environment, which broadens the spectra. This environment is fully characterized by its spectral density,</p><p>which we choose to be of the Debye form. Here, λ env is the reorganization energy of the environment that quantifies the extent to which the harmonic oscillators change equilibrium positions in going from the electronic ground to excited state, and ω c , known as the cutoff frequency, determines the timescale of relaxation of the environment. The value of ω c = 22 cm −1 was chosen to be lower than any relevant frequency characterizing the chromophore model system, whereas λ env was varied to describe both the stronger and the weaker environment coupling limits. The contribution of the environment degrees of freedom to the total response function can be exactly computed in the second-order cumulant approach,</p><p>where g env 2 (t) can be computed by substituting Eq. ( 35) into Eq. (22).</p><p>Because in our study the chromophore and environmental degrees of freedom are not directly coupled, but instead influence the energy gap independently, the total response function of the full system takes a multiplicative form, χ tot (t) = χ chr (t)χ env (t), (37)   where the exact form of the chromophore response function χ chr (t) depends on both the model system chosen for the solute degrees of freedom and the method used for evaluating the linear response function. This model for the condensed-phase environment decouples the degrees of freedom of the model system describing the chromophore from the degrees of freedom of the environment. This simplification is necessary to ensure that the spectra can be straightforwardly computed, but means that our model systems cannot capture explicit coupling of the chromophore to the condensed phase environment. Given that the improved treatment of these types of strongly coupled systems is one of the main strengths of the ensemble, E-ZTFC, and cumulant approaches over the FC approach, we note that the model systems studied in this work can be considered the best case scenario for the FC approach as applied to condensed phase systems.</p><p>Although in our models the chromophore and the condensed-phase environment are decoupled, we emphasize that the distinction between chromophore and environment in our work need not align with the solute and solvent. Indeed, from a more general perspective, many of the degrees of freedom that couple to the energy gap operator in the condensed phase, whether arising from the solute or solvent, may be captured in the linear response limit by an effective BOM. In contrast, there may be degrees of freedom arising either from the solute (e.g., intramolecular solute vibrations and torsions) or solute and solvent (e.g., hydrogen bond vibrations) whose coupling to the energy gap operator may be of a different functional form. Because this linear response description is implicitly captured by methods such as the ensemble, E-ZTFC, and cumulant approaches, our models provide rigorous and representative tests of these methods. In contrast, methods such as the FC approach, which require the selection of degrees of freedom for which one can perform a well defined normal mode analysis, explicitly perform expansions of the PESs and therefore do not necessarily include the linear response limit of realistic models. Hence, the models chosen here provide tests that only probe the best case scenario for the FC approach as applied to condensed phase systems.</p><!><p>In this work, we focus on three different model chromophore systems to analyze the effect of different approximations on the computed absorption spectrum. We first focus on the GBOM as a simple example of a system where the energy gap fluctuations do not follow Gaussian statistics. We then explore the effects of anharmonicity in the PES by studying a 1D Morse oscillator. Finally, we focus on the phenolate ion in vacuum and compute its absorption spectrum for the S 0 → S 1 transition.</p><!><p>Figure 2 shows the computed spectra for the two-mode GBOM for which the FC method yields the exact solution. In the weaker coupling limit (Fig. 2a), the ensemble approach fails to capture any of the vibronic fine structure and the pronounced asymmetry; it is also considerably too narrow. The lack of vibronic fine structure is due to the neglect of any dynamic contributions to the spectrum, whereas the lack of spectral width arises from neglecting the zero-point energy motion of the nuclei with classical sampling of configurations. 22,47 By including some dynamical effects the E-ZTFC approach correctly captures the asymmetry of the spectrum, but produces an overly broad spectrum. This feature of E-ZTFC has previously been shown to arise from the double-counting of nuclear degrees of freedom, which is most significant in the case of weaker coupling to the environment. 22,47 The second-order cumulant approach correctly captures the asymmetry of the spectrum and the vibronic fine structure. However, the peaks are red-shifted and broadened compared to the exact solution. The errors in the second- order cumulant approach result from the non-linear contribution to the fluctuations of the energy gap operator arising from both the Duschinksy rotation and the mismatch between ground and excited state frequencies. For this system, including third-order contributions to the cumulant expansion brings the spectrum into essentially quantitative agreement with the exact result.</p><p>As the coupling to the environment is increased (see Fig. 2 b), the vibronic fine structure is subsumed by contributions from the environment, leading to a single broad feature that is captured by all methods. In this limit all methods yield reasonable agreement with the exact result although the ensemble spectrum still produces a spectral shape that lacks asymmetry and the E-ZTFC approach slightly over-broadens the spectrum. The better performance of the second-order cumulant approach in the stronger chromophore-environment coupling regime is because the environment is composed of BOMs, for which the second-order cumulant method is exact. Hence, as the chromophore-environment interaction is increased, there is a larger contribution to the total spectrum from the environment in comparison to the solute, and the deviation of the second-order cumulant from the exact result decreases.</p><p>Given the strong performance of the cumulant approaches and their comparative lack of use in previous calculations of optical absorption spectra, it is worth assessing where they might work well and where they might fail in the context of the GBOM. The second-and thirdorder cumulant approaches rely on two major approximations: (i.) that the cumulant expansion can be truncated at a low-order and (ii.) that the quantum correlation functions, Eqs. ( 24) and (26), can be calculated or accurately approximated. The GBOM provides an apt test for both of these approximations because its exact spectrum can be calculated and because its harmonic nature allows for the exact calculation of the quantum correlation functions necessary to construct the cumulant spectra.</p><p>Figure 3 shows the effect of either increasing the size of the Duschinsky rotation or the difference between the ground and excited state frequencies relative to the parameters in Fig. 2. In both cases the performance of both the second-and third-order cumulant approaches degrade significantly, even leading to negative contributions to the spectrum when using the latter approach. This breakdown of the second-and third-order cumulant approaches is most marked in the weaker environment coupling regime (panels a and b of Fig. 2), although these methods still capture some of the vibronic fine structure that remains completely absent in the ensemble and E-ZTFC methods. Similar unphysical behavior in higher order cumulant approaches has been previously observed in the context of linear spectra computed for 1D anharmonic potentials. 99 Results presented here demonstrate that these breakdowns can occur even for fully harmonic PESs, as long as Duschinksy rotations and changes in frequencies between ground and excited state PES introduce sufficiently large nonlinear terms in the energy gap fluctuations. The failures of the cumulant approaches become less apparent in the stronger environment coupling regime. The reason for this, as explained above, is that the contribution arising from the environment, which is modeled using an environment composed of BOMs for which the second-order cumulant captures the exact contribution to the spectrum, dominates the deviant contribution due to the GBOM. This improved agreement between the spectra obtained using the second-and thirdorder cumulants and the exact result is also reflected in the increasingly Gaussian statistics of the energy gap fluctuation operator, shown in Fig. 4.</p><p>The shift toward increasingly Gaussian statistics in going from the weaker to stronger coupling regimes is evident in the reduction of the skewness parameter of the energy gap fluctuation operator, γ qm and γ cl , in Fig. 4 panels a. to c. and b. to d., of both the classical and quantum distributions. The skewness parameter, which is zero for a system with truly Gaussian statistics, is a way of characterizing non-Gaussian behavior. It is defined as the third moment of the fluctuation operator,</p><p>where, as defined in Sec. II, the angular brackets, ... , denote averaging with respect to the thermal distribution of the nuclear configurations of the ground state PES and σ 2 δU = δU 2 is the variance of the energy gap fluctuation. We note that the skewness parameter of the quantum distributions is higher due to the zero-point energy, which allows for fluctuations further away from the equilibrium nuclear configuration, where the non-linear component of the energy gap operator becomes more important. As we demonstrate below, the size of γ can be used as a useful diagnostic for how well the second-and third-order cumulant approaches are likely to perform for a realistic system. The origins of the successes and failures of the cumulant approaches can also be assessed by considering their approximations to the lineshape function, G n (t), used to construct the absorption spectra from Eq. ( 19). These can be compared to the exact lineshape function, G ∞ (t), by calculating the exact response function, χ, and then inverting Eq. (19). The resulting expression for the exact lineshape function for the GBOM is given in SI Eq. (S20). Because the total response function of the chromophore and environment given in Eq. ( 37) takes a product form, the contributions from the chromophore and the environment to the total lineshape function are additive,</p><p>Energy (eV) where G chr ∞ (t) and G env ∞ (t) denote the contributions to the exact lineshape function arising from the chromophore and the environment, respectively.</p><p>Figures 5 a. and b. show the exact Re[G chr ∞ (t)], the second-order cumulant approximation Re[G chr 2c (t)], and the third-order cumulant approximation Re[G chr 3c (t)] using quantum correlation functions for the two model system parameterizations used in Fig. 3.</p><p>Because Re[G chr (t)] arises from just a few modes of the GBOM representing the chromophore, it results in oscillatory behavior. The second-and third-order approximations are accurate for short times, but degrade at long times. In particular, the third-order cumulant improves agreement with the exact result over second order at short times, but becomes more ill-behaved at longer times. However, coupling to a condensed phase environment can reduce the impact of deviations of the approximate lineshape function for the chromophore at long times on the resulting absorption spectrum. For example, in the present work the environment is represented by a broad distribution of BOMs, resulting in Re[G env ∞ (t)] being a monotonically increasing function, and therefore creating an exponentially decreasing e −Re[G env (t)] (green lines in Fig. 5). The rate of this decay increases as the coupling of the energy gap to the environment increases, leading to the faster decay of the solid green line (∼25fs) in comparison to the light green line (∼75fs) corresponding to stronger and weaker environmental coupling, respectively. This exponentially decreasing environmental component of the response function, e −G env (t) , forces the total response function in Eq. (39) to zero thus setting an upper limit to the timescale for which the second-and third-order cumulant Real part of the lineshape function for the twochromophore model system as computed in the second-order cumulant (red) and third-order cumulant (blue) approaches in comparison with the exact solution (black). The real part of the total response function due to environment degrees of freedom is shown in green (dark and light green correspond to the stronger and weaker coupling regimes, respectively). The top panels correspond to the cumulant approach based on the exact quantum correlation functions, whereas the bottom panel shows lineshape functions constructed from purely classical correlation functions.</p><p>approximations to G chr ∞ (t) are required to be accurate to capture the spectrum correctly. Although stronger coupling to the environment can diminish the importance of deviations from the exact lineshape, even for a system as simple as the GBOM, the third-order cumulant can lead to instabilities. One example is the case where there is no mode mixing and the frequencies of the excited state PES are larger than those on the ground state PES, where, as shown in SI Sec. II D 3, the third-order cumulant can give unstable results, in contrast to the second-order cumulant, which is always stable under those conditions.</p><p>Although quantum correlation functions are easily calculable for harmonic models, such as the GBOM, calculating the exact quantum dynamics of general systems remains an significant challenge. Hence, it is necessary to approximate the correlation functions of the energy gap fluctuation necessary to construct the second-and third-order lineshape functions, G 2c (t) and G 3c (t). As discussed in Sec. II B 3, the conventional approximation for the second-order cumulant approach is to construct the requisite single-time quantum correlation functions from their classical analogues multiplied in Fourier space by a quantum correction factor. To obtain the multitime correlation functions necessary to construct g 3 (t) for the third-order cumulant approach, we employed the recently derived 85   third-order approximations to the lineshape function. At short times, the approximate G chr 2c (t) and G chr 3c (t) agree well with their quantum analogues in 5a. and b., and, as such, also agree with the exact result. However, the approximate G chr 2c (t) and G chr 3c (t), particularly the latter, are significantly more ill-behaved at longer times. G chr 2c (t) is captured more accurately because the quantum correction factor used here is exact for single-time linear correlation functions of the positions of a harmonic Hamiltonian and these types of correlation functions dominate g 2 (t) for the GBOM (See SI Sec. II D). Calculation of the additional g 3 (t) contribution needed for the third-order cumulant requires nonlinear correlation functions of the coordinates, which are more challenging to approximate using classical input. For the GBOM, the approximate quantum correction factor introduces a divergent term in G chr 3c (t), which cancels exactly when using the exact quantum correlation function for any choice of model system parameters (See SI Sec. II D 3), but gives rise to the unphysical behavior of G chr 3c (t) at long timescales shown in Figs. 5 c. and d.</p><p>Despite the differences in the quantum and approximate correlation functions at longer times, to obtain an accurate spectrum the time over which the correlation function needs to be accurately captured is limited to short timescales because of the coupling to the environment. Over these timescales the discrepancies between the lineshape function constructed from the exact quantum correlation function and the classical correlation function with quantum correction factor are small. Section III of the SI demonstrates that these arguments also apply to the imaginary parts of the second-and third-order cumulant approximations to the lineshape functions. The deviations from the exact quantum cor-relation functions manifest as minor discrepancies in the spectra generated from the second-and third-order cumulant approaches using exact quantum input, as Fig. 6 shows. Given that the exact quantum correlation functions are generally only accessible for a few simplified model systems, this is a reassuring result for applying the cumulant approach to more realistic systems. In the following sections, we employ the cumulant approach as constructed from purely classical correlation functions, assessing its performance for both the anharmonic Morse oscillator and the phenolate ion.</p><!><p>We now turn our attention to the Morse oscillator model that contains anharmonicity in both the ground and excited state PESs. Unlike the GBOM model, for which the harmonic FC method provides the exact spectrum, the Morse oscillator provides a stringent test for methods that employ harmonic approximations, such as the FC and E-ZTFC schemes. Although the second-order cumulant approach does not explicitly require harmonic PESs, the fact that it is exact for systems where the statistics of the energy gap fluctuation are Gaussian, such as the BOM, raises fundamental questions about its ability to capture anharmonic effects. Figure 7 shows the spectra calculated for the Morse oscillator model parameterized to the carbon monoxide molecule, CO, immersed in a BOM environment coupled weakly to the energy gap (see SI Sec. IV). As was shown Sec. V A, the case of weaker coupling to the environment is a more challenging test of these methods, so here we examine weaker coupling to the environment and provide results for stronger coupling to the environment in the SI (Sec. IV). Here, we used classical sampling for the ensemble and E-ZTFC methods and classical dynamics to generate the correlation functions necessary for the second-and third-order cumulant approaches.</p><p>Figure 7 a. shows that the FC method relying on the harmonic expansion of the ground and excited state PESs accurately captures the vibronic progression at lower energies, but fails to capture the high energy peaks. The FC method maps the Morse oscillator model to an effective GBOM, for which it then calculates the exact optical spectrum. However, the sharp increase in the harmonic potentials in the GBOM overestimate the increase in energy with the extension of the CO bond and hence neglect the large extensions allowed by the Morse oscillator that give rise to the high energy peaks. Consistent with the spectrum for the GBOM, the ensemble method misses the asymmetric vibronic progression. The E-ZTFC scheme captures the width and asymmetry of the exact spectrum more accurately than the ensemble method, but, because it relies on the zero-temperature limit of the harmonic approximation to the FC method, also fails to capture the higher energy features. Figure 7 b. considers a more strongly anharmonic PES, reveal-ing that the FC method is insensitive to the increase in anharmonicity because the GBOM to which the Morse oscillator system is mapped is identical to the less anharmonic Morse oscillator. The ensemble and E-ZTFC methods result in broader peaks, albeit not sufficiently broad to capture the breadth of the exact spectrum.  d. show the second-and third-order cumulant spectra for the Morse oscillator parameterization of the CO molecule and its more anharmonic variant, respectively. In contrast to the FC spectra, the secondand third-order order cumulant approaches capture the higher energy vibronic progression of the exact spectrum. Both the second-and third-order order cumulant schemes overestimate the separation between subsequent vibronic peaks, with the third-order cumulant providing a closer estimate. The skewness parameter, γ α CO = −0.104 and γ 1.8α CO = +0.071, for the energy gap fluctuations in both the CO Morse system and its more anharmonic variant are small. This small value indicates that the second-and third-order cumulant approaches are likely to perform well, reinforcing the observations made in Figs. 7 c. and d. As we show in Fig. S3 (SI Sec. IV), when the coupling to the environment is increased, the performance of the ensemble and E-ZTFC schemes improves, although they still perform much worse than the second-and thirdorder cumulant approaches which match the exact spectra very closely. In contrast, the FC results closely resemble those of E-ZTFC, missing the pronounced asymmetry that the second-and third-order cumulant approaches capture.</p><p>Because the second-order cumulant approach is exact for systems with energy gap fluctuations displaying Gaussian statistics, such as the BOM, it may seem that the method could potentially only capture harmonic effects. However, the results for the Morse oscillator demonstrate that the second-order cumulant captures features of the spectrum beyond those present in the FC method, which corresponds to the harmonic expansion of the PESs. Indeed, the second-order cumulant maps the fully anharmonic Morse oscillator system to an effective BOM characterized by harmonic potentials with their minima separated by a smaller distance than that of the original Morse system. This approximate mapping allows the second-order cumulant to mimic the spectral features arising from anharmonic effects in the exact spectrum. It also explains why, for the cumulant approach, the onset of the first vibronic peak is at higher energy in Fig 7  c: although in this case the cumulant spectrum correctly captures the average energy gap, thereby fixing the center of the spectrum, the spectral shape is that of the effective BOM with a reduced shift vector, resulting in the 0-0 transition being shifted up in energy by approximately 0.17 eV compared to the exact and the FC 0-0 transition (see SI Sec. IV A for a more detailed analysis of the spectral shift and BOM mapping).  To mimic the broadening effects of the condensed phase, the phenolate ion was coupled to a BOM-type environment characterized by λenv = 0.073eV. Spectra are shifted such that the vertical excitation energy computed at the minimum of the ground state PES is located at ω = 0 eV.</p><!><p>We now consider the spectrum arising from the S 0 → S 1 transition for a fully atomistic representation of the phenolate ion. In this case, as detailed in Sec. III D, we performed ab initio molecular dynamics at the DFT level to thermally sample the ground state configurations of the chromophore required for the ensemble and E-ZTFC methods and to provide the correlation functions needed for the cumulant approaches, which was then immersed in a weakly coupled BOM environment. As in previous examples, Fig. 8 shows that the ensemble and E-ZTFC approaches yield spectra devoid of the vibronic progression. The FC, second, and third-order cumulant schemes all give similar results. The small skewness parameter of the energy gap fluctuations, calculated using classical sampling on the ab initio PES, suggest that the second-and third-order cumulant schemes can be expected to provide a close approximation to the exact spectrum. Further, the energy gap fluctuations of the effective GBOM arising from the normal mode analysis of the phenolate PESs also display nearly Gaussian statistics, with skewness parameters of γ = 0.047 and γ = 0.088 for the classical and quantum distributions, respectively (see Fig. 9). The close agreement of the FC results with those of the second-and third-order cumulant approaches suggests both that all three methods likely provide an accurate approximation to the exact spectrum and also that the harmonic potential approximation underlying the FC method is a good approximation for this system. However, despite the close agreement between the optical spectra obtained using these methods, the second-and third-order cumulant approaches likely capture subtle anharmonicities that the FC method misses. This can be observed in the differences between the spectral density obtained from directly sampling the phenolate PESs and the spectral density obtained from its parameterization to an effective GBOM arising from normal mode analysis (see SI Sec. V).</p><!><p>Here, we have studied methods for computing linear absorption spectra that range in their abilities to address how anharmonicity, the condensed phase environment, and dynamical effects manifest in the optical absorption spectrum. In the presence of stronger environment coupling, all of the methods considered here, save the ensemble approach, are able to capture the broad band of the optical spectra when the surfaces are harmonic. When the surfaces are anharmonic, only cumulant-based approaches are able to capture the spectra. In contrast, when the environment couples weakly to the energy gap, vibronic progressions become evident and the spectral signatures of the chromophore become more prominent. In such cases, the dynamical effects that are completely neglected in the ensemble and partially ignored in the E-ZTFC methods lead to their inability to capture vibronic progressions. Instead, the FC, second-, and thirdorder order cumulant approaches capture the vibronic progressions, leading to generally qualitative and sometimes quantitive agreement with the exact spectra. We also assessed the ability to approximate the quantum correlation functions required for the second-and thirdorder cumulant approaches using classical dynamics with quantum correction factors. Due to the environmental coupling in the condensed phase, the spectrum is only sensitive to the short time (∼10s of femtoseconds) part of these correlation functions, supporting the use of such quantum correction schemes.</p><p>While the model systems studied here present a controlled set of challenges to these methods, realistic systems introduce distinct complications for which these methods may fail. For example, in a real, condensed phase system, the distinction between chromophore and environment becomes less well-defined. This presents a particular challenge for the FC approach as it requires optimized geometries of a chromophore subsystem, i.e., a part of the chromophore, the entire chromophore, or the chromophore and strongly hydrogen bonded solvent, on the ground and excited state PESs which are needed to calculate the Hessians. Therefore, calculation of the necessary Hessians can be ill-defined in the condensed phase.</p><p>Although the second-and third-order cumulant approaches have not been heavily used to calculate optical spectra in the condensed phase due to their high computational cost, they provide distinct advantages: (i) unlike the ensemble and E-ZTFC methods, they capture dynamical effects such as vibronic progressions, (ii) they can easily include the effect of the environment on the energy gap, circumventing the problems in applying the FC method to condensed phase systems, and (iii) they are applicable to anharmonic systems for which the statistics of the energy gap fluctuation are (nearly) Gaussian. As we have shown, using the same data required to calculate the ensemble or E-ZTFC spectra, one can calculate the skewness parameter, γ of the energy gap fluctuations, which serves as a diagnostic of how well the second-and third-order cumulant approaches will perform. Indeed, as we have demonstrated, when γ < 0.2 the secondand third-order cumulant approaches can provide close agreement with the exact spectrum. With the advent of advances in excited state electronic structure theory algorithms and high-performance computing capabilities, the use of cumulant methods, which require many thousands of excited state electronic structure calculations, is now possible for condensed phase systems.</p>