Optical spectra in the condensed phase: Capturing anharmonic and vibronic features using dynamic and static approaches

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–5 of a material constitutes an essential first step in understanding its nonlinear spectroscopy.^6–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–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 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.

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–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).¹ The choice of these distributions is motivated by the inhomogeneous (Gaussian) and homogeneous (Lorentzian) limits corresponding to slowly and rapidly evolving nuclear coordinates, respectively,²³ coupled weakly to a classical vertical energy gap obeying Gaussian statistics.²³ The ensemble method^{13,16,24–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 time scales.⁶ It is important to note that although molecular dynamics (MD) 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 that leads to features such as vibronic fine structure.

The Franck-Condon (FC) method^28–40 and related schemes^{2,3,5,41–43} 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. Although recent applications of the FC method can approximately account for anharmonicity in the ground and excited state PESs,^42,44 these applications rely on computationally expensive third derivatives of the energy with respect to nuclear positions, making their broad application challenging. As such, the common implementation of the FC approach employs the harmonic approximation to the ground and excited state PESs to estimate the eigenstates and therefore poorly describes 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.⁴⁵ 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 invalidate 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.^1,22,46 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.

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,47 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 the anharmonic effects and the dynamical effects that lead to vibronic spectral 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 require only the calculation of the equilibrium energy gap fluctuation (single time) autocorrelation function.^6,48 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. 47 and 49). The second-order cumulant scheme has been heavily exploited to reproduce and elucidate the vibrational spectroscopy of probes in the condensed phase^9,50–52 and has also recently found its use in directly simulating electronic spectra in the condensed phase.^53–55 Beyond linear spectroscopy, it can also be used to investigate nonlinear spectroscopy and energy transfer in molecular systems.^6,48 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,⁵⁶ amide and carbonyl in water,⁵⁷ and the OH stretch in deuterated water.⁵⁸ 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.

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.

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,

where g and e denote the electronic ground and first excited state, respectively, $Hgq^$ and $Heq^$ are the nuclear Hamiltonians corresponding to the ground and excited electronic states, respectively, and $Vegq^=Vgeq^†$ are the transition dipole operators, which are generally dependent on the nuclear configuration.

For this system, the linear absorption spectrum can be written as^6,48

where α(ω) is an ω-dependent prefactor whose form is dependent on the particular definition of the absorbance σ(ω) used.^6,42 For the purpose of this work, we use α(ω) = 1 such that the reported absorption spectra are directly proportional to the Fourier transform of χ(t). The linear response function χ(t), in atomic units, is given by

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-ordered exponentiated integral of the energy gap operator,

where $Vg(q^)$ and $Ve(q^)$ are the ground and excited state PESs, respectively, and $ωeg0$ 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)=eiĤg(q^,p^)tO(p^,q^)e−iĤg(q^,p^)t$⁠.

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,

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 Secs. II A and II B, 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.

The approaches presented in this section rely on the assumption that the nuclear motions are slow in comparison to the time scale 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:

where

is the representation of the canonical nuclear operator in real space and R denotes the choice of coordinate system used to describe the nuclear positions. 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:

The second line of Eq. (6) represents a convenient starting point for a range of exact to approximate treatments. For example, by using path integral molecular dynamics (PIMD)⁵⁹ 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.

The simplest static approach can be obtained by setting $f(R)=δR−Rg0$ in Eq. (7), which is exact in the classical limit at zero temperature, and where $Rg0$ corresponds to the energy minimum on the ground state PES. Upon substitution into Eq. (8), this yields

Thus, the absorption spectrum reduces to a δ-function placed at a vertical excitation energy computed for a single representative nuclear configuration $Rg0$⁠. To generate a spectrum with a finite width, $δ(ω−U(Rg0))$ is often replaced by a Gaussian function, where the standard deviation is chosen in a phenomenological way.¹ Because finding a single well-defined global minimum of the ground state PES is impractical in disordered condensed phase systems, $Rg0$ is often computed for a chromophore embedded in continuum solvent.^4,20,60 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 continuum solvation model. However, the approach requires only a single evaluation of the vertical energy gap, $U(Rg0)$⁠, making it the only tractable method if high-level wavefunction schemes are required to compute reliable excited state energies for the system of interest.

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⁵⁹ or by sampling from the Wigner distribution.⁶¹ 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)=∫dPe−H(P,R)Z$ which corresponds to the classical configurational distribution function.

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 lead to an inconsistency between the vertical energy gap evaluation and the underlying sampling.^53,62,63 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¹ 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.

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.

The FC method^28–40 generally 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,^64,65

Several approaches exist that relax the Condon approximation by allowing V_eg to vary linearly with nuclear degrees of freedom,^5,41–43 but these Herzberg-Teller type 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

Here, {|g⟩ ⊗ |ν_g⟩, |e⟩ ⊗ |ν_e⟩} are the vibronic eigenbases of the Hamiltonian, Eq. (1a), within the harmonic approximation such that $Hg(q)|νg=Eνg|νg$ and $He(q)|νe=Eνe|νe$⁠, $Ωνe,νg=Eνe−Eνg$ is the energy gap between the ν_gth and ν_eth vibrational 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

As in the case of the vertical excitation gap and ensemble methods, the delta functions in Eq. (12) are often broadened using Gaussian distributions.

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 vibrational state of the ground PES, i.e., $ρνg→δ0g,νe$ which leads to the following expression for the absorption spectrum:

and its associated zero-temperature FC (ZTFC) response function,

The FC approach has been shown to yield accurate absorption spectra in a variety of systems,^2,66–68 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.⁶⁹ 

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.²² The approach is aimed at retaining the main advantages of the ensemble approach, such as the explicit sampling of chromophore-environment 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.

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 ensemble method, Eq. (8). The second component, consisting of the ZTFC response, Eq. (14), is temperature-independent and accounts for the vibronic contribution to the spectrum in the zero temperature limit, i.e., those arising from the ground vibrational state. For many high-frequency intramolecular vibrations, only the ground vibrational state is populated and the zero temperature limit agrees with the room temperature result. Hence, the E-ZTFC response function can be constructed as follows:

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 time scale 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,

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,

where $Rg0$ is the nuclear configuration at the minimum of the ground state PES and $U(Rg0)$ is the vertical excitation energy at this position.

To obtain the quantum limit of Eq. (16), i.e., when using PIMD to sample the ground state PES in the ensemble method, we exploit the fact that within the FC spectrum, the most general form of the parameterized 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 supplementary material. Thus, the quantum limit of Eq. (16) corresponds to the β → ∞ and t → 0 limits of the GBOM response function, Eq. (S17) in Sec. II B of the supplementary material, for the parameterized version of the GBOM arising from the normal mode analysis used in the FC approach.

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 $χZTFCinh(t)$ provided in Sec. II C of the supplementary material. Substituting Eqs. (15) and (16) into Eq. (2), we obtain

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.

The cumulant scheme^6,47 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)−ωegav$⁠,

where

is the mth order approximation to the full line shape function G_∞(t). Here, g_n(t) is the nth order cumulant and $ωegav=⟨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,48

The mth order cumulant approximation to the response function can then be expressed in terms of the mth order approximation to the line shape function,

For the second-order cumulant, i.e., g₂(t) = G₂(t), where the expansion is truncated at second-order, the line shape function takes the form

The second-order line shape function depends on the spectral density, $J(ω)$⁠, which can be expressed as follows:

where θ(ω) is the Heaviside function and

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, such as the GBOM.

Although many schemes have been developed to approximate equilibrium time correlation functions, the standard approximation in the application of the second-order cumulant approach to general anharmonic systems is to use the classical version of C_δU(t) to construct an approximate quantum analog.^6,70–75 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 Kubo-transformed analog.^76,77 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δUkubo(t)$⁠, leading to the following expression for the spectral density:

The specific choice of quantum correction factor employed in Eq. (25), while not unique,⁷¹ has been widely used in the context of computing linear spectra of chromophores in complex environments^53–55 and has shown good performance in constructing spectral densities for studying the energy transfer in multichromophore pigment-protein systems.^78–82 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δUcl(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 line shape 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).

The second-order cumulant approximation therefore relies on two central approximations. First is the second-order 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 second-order 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.

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, $g3(t)=g3CδU(3)(τ1,τ2)(t)$⁠, 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 Sec. I of the supplementary material.

A major difficulty that the third-order cumulant presents is the need to calculate a multitime quantum correlation function,

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 multitime quantum correlation functions in terms of their classical counterparts. Recently, an approximate link has been derived⁸³ to connect the Fourier transform of the double Kubo-transformed quantum correlation function to its classical counterpart,

where $ω¯=ω+ω′$⁠. Here, we assess this recent development in the context of linear optical spectroscopy.

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 in the E-ZTFC approach due to decomposing the spectrum into temperature-independent dynamic and finite-temperature 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 time step 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 continuum model.

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 of 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 Sec. IV is for a BOM-type environment.

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,

As defined at the beginning of Sec. II, $ωeg0$ is the energy difference between the minima of the multidimensional ground and excited state PESs or adiabatic energy gap, $p^j$ and $q^j$ are the mass-weighted momentum and coordinate operators (i.e., $p^j=P^j/mj$ and $q^j=mjQ^j$⁠), ω_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 line shape function. Here, the energy gap fluctuation operator is linear in the oscillator coordinates, $δU(q)=U(q)−ωegav=∑jωj2Kjq^j$⁠, where $ωegav=ωeg(0)+λ$ with $λ=12∑jωj2Kj2$ is called the reorganization energy.

Harmonic oscillators are sufficiently simple to allow analytical solutions to their quantum dynamics; the exact expression for the spectral density, $J(ω)$⁠, in Eq. (23)⁶ is

Because the second-order cumulant expansion is exact for the nuclear Hamiltonian in Eqs. (28a) and (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.

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,

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,⁸⁴ 

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,^42,43,85,86 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 Sec. II of the supplementary material). In particular, the harmonic form of the PES means that exact expressions for the one- and two-time quantum correlation functions can be derived (Sec. II D of the supplementary material), 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.^54,55

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 nonlinearities 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 states 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

and a value of ϕ = 10° is used for the GBOM results presented in this work unless otherwise noted. This choice is again informed by typical parameters for realistic systems, where some mode mixing occurs, but the Duschinsky matrix is found to be diagonally dominant.⁴² 

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

where $D=k2α2$ is the well depth, k is the force constant at the minimum of the well, α is the anharmonicity parameter (with the limit of α → 0 recovering the harmonic oscillator upon Taylor expansion of the exponential), 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 surfaces to a scan over the bond stretching mode of the carbon monoxide molecule. The calculations were performed in Gaussian⁸⁷ using density functional theory (DFT) to describe the ground state and time-dependent DFT (TDDFT) in the Tamm-Dancoff approximation⁸⁸ to describe S₁ excited state PES with the CAM-B3LYP⁸⁹ functional and the 6-31+G^*90 basis set. Further details of the fit and the resulting Morse potential parameters can be found in Sec. V of the supplementary material. 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², 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,⁹¹ 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 the 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.

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

where μ is the effective mass of the vibrational mode. To compute the spectra for this system in the cumulant approach, we performed classical MD at a temperature of 300 K on the one-dimensional ground state PES using the i-PI code.^92,93 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 Sec. V of the supplementary material.

Finally, we move away from model systems to study the S₀ → S₁ transition of the phenolate (C₆H₅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 line shape. 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,⁸⁷ where the Tamm-Dancoff approximation⁸⁸ is applied to the excited state TDDFT calculations. The shift vector K and the Duschinsky 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 (AIMD) simulation of phenolate in vacuum at a temperature of 300 K using the TeraChem code.^94,95 The same basis set and DFT functional as for the normal mode calculations were used to perform the MD. Detailed simulation parameters for the AIMD can be found in Sec. VI of the supplementary material. Three independent 50 ps trajectories were generated from the simulations. From these trajectories, the chromophore coordinates were extracted every 2 fs and the S₁ vertical excitation energy was calculated using the TeraChem code⁹⁶ 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.

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,

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 the excited state, and ω_c, known as the cutoff frequency, determines the time scale of relaxation of the environment. The value of ω_c = 22 cm⁻¹ 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,

where $g2env(t)$ can be computed by substituting Eq. (35) into Eq. (22).

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,

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.

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.

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₀ → S₁ transition. All systems are simulated at a temperature of 300 K. For the GBOM, cumulant spectra are evaluated using analytical expressions for the exact one- and two-time quantum correlation functions, unless specified otherwise. For the Morse oscillator and the phenolate ion, all cumulant spectra are evaluated by approximately reconstructing the quantum correlation functions from their classical counterparts obtained directly from the MD, using Eq. (25) for the second order and Eq. (27) for the third order cumulant contribution. Ensemble and E-ZTFC spectra are computed using classical sampling on the ground state PES, obtained either from the relevant analytical expressions for the response function of the GBOM derived in Sec. II of the supplementary material or through classical MD sampling of nuclear configurations for the Morse oscillator and the phenolate ion.

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. 2(a)], 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,46 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,46 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 nonlinear contribution to the fluctuations of the energy gap operator arising from both the Duschinsky 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.

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.

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 third-order 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 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.

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 [Figs. 2(a) and 2(b)], 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.⁹⁷ Results presented here demonstrate that these breakdowns can occur even for fully harmonic PESs, as long as Duschinsky 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 third-order cumulants and the exact result is also reflected in the increasingly Gaussian statistics of the energy gap fluctuation operator, shown in Fig. 4.

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 Figs. 4(a)–4(c) and Figs. 4(b)–4(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,

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 $σδU2=⟨δU2⟩$ 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 nonlinear 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 line shape function, G_n(t), used to construct the absorption spectra from Eq. (19). These can be compared to the exact line shape function, G_∞(t), by calculating the exact response function, χ, and then inverting Eq. (19). The resulting expression for the exact line shape function for the GBOM is given in Eq. (S20) of the supplementary material. 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 line shape function are additive,

where $G∞chr(t)$ and $G∞env(t)$ denote the contributions to the exact line shape function arising from the chromophore and the environment, respectively.

Figures 5(a) and 5(b) show the exact $Re[G∞chr(t)]$⁠, the second-order cumulant approximation $Re[G2cchr(t)]$⁠, and the third- order cumulant approximation $Re[G3cchr(t)]$ using quantum correlation functions for the two model system parameterizations used in Fig. 3. 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 line shape 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[Genv(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 dark green line (∼25 fs) in comparison to the light green line (∼75 fs) corresponding to stronger and weaker environmental coupling, respectively. This exponentially decreasing environmental component of the response function, $e−Genv(t)$⁠, forces the total response function in Eq. (39) to zero thus setting an upper limit to the time scale for which the second- and third-order cumulant 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 line shape, 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 Sec. II D 3 of the supplementary material, the third-order cumulant can give unstable results, in contrast to the second-order cumulant, which is always stable under those conditions.

Although quantum correlation functions are easily calculable for harmonic models, such as the GBOM, calculating the exact quantum dynamics of general systems remains a significant challenge. Hence, it is necessary to approximate the correlation functions of the energy gap fluctuation necessary to construct the second- and third-order line shape 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 analogs multiplied in Fourier space by a quantum correction factor. To obtain the multitime correlation functions necessary to construct g₃(t) for the third-order cumulant approach, we employed the recently derived⁸³ quantum correction factor connecting quantum and classical double-Kubo transformed correlation functions. Figures 5(c) and 5(d) show the second- and third-order approximations to the line shape function. At short times, the approximate $G2cchr(t)$ and $G3cchr(t)$ agree well with their quantum analogs in Figs. 5(a) and 5(b), and, as such, also agree with the exact result. However, the approximate $G2cchr(t)$ and $G3cchr(t)$⁠, particularly the latter, are significantly more ill-behaved at longer times. $G2cchr(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₂(t) for the GBOM (see Sec. II D of the supplementary material). Calculation of the additional g₃(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 $G3cchr(t)$⁠, which cancels exactly when using the exact quantum correlation function for any choice of model system parameters (see Sec. II D 3 of the supplementary material) but gives rise to the unphysical behavior of $G3cchr(t)$ at long time scales shown in Figs. 5(c) and 5(d).

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 time scales because of the coupling to the environment. Over these time scales, the discrepancies between the line shape function constructed from the exact quantum correlation function and the classical correlation function with quantum correction factor are small. Section IV of the supplementary material demonstrates that these arguments also apply to the imaginary parts of the second- and third-order cumulant approximations to the line shape functions. The deviations from the exact quantum correlation functions manifest as minor discrepancies in the spectra generated from the second- and third-order cumulant approaches using exact quantum input, as shown in Fig. 6. 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 Secs. V B and V C, 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.

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 Sec. V of the supplementary material). As shown in 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 Sec. V of the supplementary material. 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.

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 overestimates the increase in energy with the extension of the CO bond and hence neglects the large extensions allowed by the Morse oscillator that gives rise to the high energy peaks. Consistent with the ensemble 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, revealing 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.

Figures 7(c) and 7(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 second- and 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 7(d). As shown in Fig. S3 (Sec. V of the supplementary material), 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 third-order cumulant approaches that 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.

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 Sec. V A of the supplementary material for a more detailed analysis of the spectral shift and BOM mapping).

We now consider the spectrum arising from the S₀ → S₁ 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, suggests that the second- and third-order cumulant schemes can be expected to provide a close approximation to the exact spectrum. Furthermore, 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 Sec. VI of the supplementary material).

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 third-order order cumulant approaches capture the vibronic progressions, leading to generally qualitative and sometimes quantitative agreement with the exact spectra. We also assessed the ability to approximate the quantum correlation functions required for the second- and third-order 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 (∼10 s of femtoseconds) part of these correlation functions, supporting the use of such quantum correction schemes.

Although 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.

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 as follows: (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 second- and 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.

See supplementary material for expressions for the second- and third-order cumulants, GBOM response functions for the various spectroscopic methods, numerical parameters used for simulating the two-mode GBOM spectra, graphs detailing the imaginary component of the GBOM line shape function, and additional information about the 1D Morse oscillator and phenolate simulations.

This work was supported by the Department of Energy, Office of Basic Energy Sciences CTC and CPIMS programs, under Award No. DE-SC0014437. This work used the XStream and MERCED computational resources supported by the NSF MRI Program (Grant Nos. ACI-1429830 and ACI-1429783).