A partially linearized spin-mapping approach for simulating nonlinear optical spectra

We present a partially linearized method based on spin mapping for computing both linear and nonlinear optical spectra. As observables are obtained from ensembles of classical trajectories, the approach can be applied to the large condensed-phase systems that undergo photosynthetic light-harvesting processes. In particular, the recently derived spin-PLDM method has been shown to exhibit superior accuracy in computing population dynamics compared to other related classical-trajectory methods. Such a method should also be ideally suited to describing the quantum coherences generated by interaction with light. We demonstrate that this is indeed the case by calculating the nonlinear optical response functions relevant for the pump--probe and 2D photon-echo spectra for a Frenkel biexciton model and the Fenna--Matthews--Olsen light-harvesting complex. One especially desirable feature of our approach is that the full spectrum can be decomposed into its constituent components associated with the various Liouville-space pathways, offering a greater insight beyond what can be directly obtained from experiment.


I. INTRODUCTION
Nonlinear optical spectroscopy is a powerful tool for elucidating the exciton dynamics of condensed-phase systems. 1 While linear spectroscopic techniques are essentially limited to determining the exciton energies and dephasing rates, nonlinear methods give access to a plethora of extra information. 2,3 For example, 2D optical photon-echo spectroscopy has previously been used to measure exciton couplings, 2,4,5 excitonic energy-transfer pathways and timescales, 2,4-7 exciton coherences 7,8 and dipole-moment orientations. 9 By tuning the polarization of each of the pulses, the relative intensity of peaks within 2D photon-echo spectra can be altered, providing a way to disentangle and assign peaks within crowded regions of the spectra, which are often present for large systems. 4,10 Nonlinear approaches are also useful for removing the effects of inhomogeneous broadening. 1 Because of these advantages, 2D optical photon-echo spectroscopy has become the experimental tool of choice for investigating the excitonic energy-transfer processes within photosynthetic light-harvesting systems, where researchers are still trying to understand their impressive ability to harvest sunlight with near perfect quantum efficiency. For theory to help interpret and give further insight to such measurements, there is a need to develop approaches for simulating the processes and calculating spectroscopic quantities which can be applied to systems with a large number of degrees of freedom.
A common way of computing such spectroscopic quantities is by using the optical response function approach. 1 Here the field-matter interaction is treated perturbatively, leading to a set of multi-time correlation functions which involve only the field-free dynamics. Each speca) Electronic mail: jonathan.mannouch@phys.chem.ethz.ch b) Electronic mail: jeremy.richardson@phys.chem.ethz.ch troscopic observable is then typically obtained by various Fourier transforms of the associated optical response functions. Alternatively the so-called explicit light-field approach [11][12][13] can also be used to compute spectroscopic properties, by performing the dynamics associated with the coupling to a classical light field through the use of a time-dependent Hamiltonian. While this alternative approach has the advantage that it can describe processes involving strong light fields, the optical response function approach is nevertheless sufficient for and provides the simplest connection to many experiments. It also has the added advantage that the full signal can be separated into contributions associated with different Liouville-space pathways, offering greater insight beyond what can be directly obtained from experiment. Many methods have been developed to compute linear and nonlinear spectra within this framework. [14][15][16][17] One of the most commonly used techniques for simulating nonlinear spectra is the Hierarchical equations of motion (HEOM) approach. [18][19][20][21] Although this method enables numerically exact quantum dynamics to be applied to large condensed-phase systems containing a harmonic bath, it is unable to treat more realistic anharmonic problems. A full quantum mechanical approach is also in some sense unnecessary for treating the problems we consider in this paper, as we expect that it should be possible to simulate most relevant biological processes without taking the quantum-mechanical nature of nuclear dynamics into account. We are therefore primarily interested in classical-trajectory techniques, which can in principle be applied to systems with any form for the nuclear potential.
When modelling spectroscopic quantities involving electronic transitions between the ground and a single excited state, there are several different classical-trajectory approximations that can be used; these differ by whether the nuclear degrees of freedom are evolved on the groundstate potential, excited-state potential or some linear combination of both. 1,22 All such approaches are exact for shifted harmonic models where the ground and excited state have the same frequencies, but give different approximations in general. From previous work, it has been found that the Wigner averaged classical limit (WACL), 23 which evolves the nuclear degrees of freedom on the time-independent arithmetic mean surface of the ground and a single excited state, generally gives rise to the most accurate results. The forward-backward initialvalue representation (FB-IVR), 24 which is related to the Herman-Kluk propagator, 25 can additionally exactly reproduce spectra for systems that involve a frequency change under photoexcitation, but this approach cannot easily be converged for high-dimensional condensedphase problems and the Herman-Kluk prefactor becomes unstable when applied to anharmonic models. 26 While WACL has been relatively successful in accurately obtaining Frank-Condon spectra for certain condensed-phase systems, it cannot describe nonadiabatic transitions between coupled exciton states and so is not an appropriate method for photosynthetic lightharvesting systems. Nonetheless, it is desirable that a nonadiabatic dynamics method reduces to WACL in the case that the excited states are uncoupled. Other features that we would want from such an approach is that it is exact in the static-nuclear limit (when inhomogeneous broadening dominates), is accurate in the hightemperature limit (where classical trajectories are expected to be valid) and can be used to compute not just single-time correlation functions, but also the multi-time correlation functions needed for obtaining nonlinear spectra.
Mapping-based approaches provide a way to accurately include nonadiabatic transitions within a classical trajectory picture. The combined dynamics of the nuclei and excitons are treated as an ensemble of classical trajectories, where the dynamics of the later is described by evolving a set of continuous variables associated with a mapping space. The computational cost of such methods therefore scales with the number of degrees of freedom in the same way as standard classical molecular dynamics, allowing them to be applied to large condensed-phase systems. Originally, most mapping-based approaches employed the so-called Meyer-Miller-Stock-Thoss (MMST) mapping, 27 where the exciton subsystem is described by the single-excitation space of a set of harmonic oscillators and the dynamics are obtained by evolving the classical phase-space variables associated with this mapping space. Recently, progress has been made in improving the accuracy of such methods through the use of a resolution of the identity, 28,29 optimized zero-point energy parameters, 30,31 the generalized master equation, 32 nonadiabatic ring-polymer molecular dynamics, 33 symmetric windowing (SQC), 34 spin-mapping, 35,36 and other alternative classical mapping models. [37][38][39][40] Some of these advancements have also already been used to obtain both linear and nonlinear optical spectra. 11,[41][42][43][44] A particularly successful method for computing dynamical observables within exciton systems is the standard partially linearized density matrix (PLDM) [45][46][47][48] approach, which uses coherent states within the MMST mapping space to describe the dynamics associated with the forward and backward exciton paths separately through the use of two independent sets of mapping variables. Additionally, a procedure for the application of this method to the computation of linear and nonlinear optical spectra has already been developed, which has been used successfully for a range of systems. 49,50 Such an approach seems especially attractive, as it fulfills many of the desirable features that we would like; the approach is exact in the static-nuclear limit, can compute both single-and multi-time correlation functions directly 49,50 and for relatively short propagation times is seen to be extremely accurate in the high-temperature limit. 46 We note that the approach also reduces to WACL in the absence of diabatic couplings if so-called 'focused' initial conditions are used. Despite the many positives, there are however certain aspects of the method that could be improved.
For standard PLDM, the initial mapping variables can be sampled using either focused initial conditions or a Gaussian sampling approach, which are both defined in Ref. 45. It is in general not obvious which sampling approach is superior. Focused conditions, which limit the initial sampling space of the mapping variables to a region corresponding to the occupation of a single initial excitonic state, have typically been used when computing nonlinear spectra for reasons of improved efficiency. 50 However it is also known that the accuracy of population dynamics can suffer dramatically by using focused initial conditions instead of a Gaussian sampling approach for the mapping variables. 45 Additionally, it has been observed that even in the high-temperature limit, standard PLDM becomes inaccurate for long propagation times even when using Gaussian sampling. 46 Such an approach could therefore lead to significant errors in the nonlinear spectra when considering long delay times between the pump and probe pulses. One potential way to alleviate some of these problems is by using spin-PLDM, 46,47 which is an extension of the standard method that instead utilizes the spin-mapping space reformulated in terms of Stratonovich-Weyl kernels. 35,36 For spin-mapping, focused initial conditions are just as rigorous in the sense that this sampling of the mapping variables leads to the product of any two electronic operators still being correctly described by the mapping, 47 while for focused conditions in the MMST mapping space this is only true for products of population operators. As a result, there is no reason to expect a loss in accuracy when focused initial conditions are used with spinmapping approaches, as has been observed, and using such initial conditions still gives similar improvements in efficiency as for standard PLDM. 35,36,46 Spin-PLDM has also been shown to exhibit greater accuracy than standard PLDM for a range of different problems and in particular it is found to be numerically exact in the hightemperature limit, even for long propagation times. 46 For these reasons, spin-PLDM is expected to offer a number of key advantages for calculating nonlinear spectra over the standard method and is the only approach of those outlined above that is able to fulfill all of the desirable features that we listed earlier.
After first reviewing how linear optical spectra can be calculated with fully linearized semiclassical (LSC) mapping methods, we introduce our approach for calculating both linear and nonlinear spectra with spin-PLDM. The improved accuracy of such an approach is then demonstrated by comparing the optical absorption, fluorescence, pump-probe and two-dimensional (2D) photonecho spectra for a range of model systems calculated with both spin-PLDM and other commonly used methods. In the literature, a particular Frenkel biexciton model is often used to benchmark approximate methods for computing nonlinear optical spectra. However, we will show that a static-nuclear picture often suffices to obtain reasonably accurate results for this system, and so in this paper we will also test the methods on more difficult parameter regimes. In addition, we compute pump-probe spectra for the seven-state Fenna-Matthews-Olsen complex (FMO) model, which as well as being relevant for studying photosynthetic light-harvesting processes, is also a much more challenging benchmark system. This is in particular because in this model, there are 21 doubleexciton states needed for computing nonlinear spectra, whereas biexciton models have only one. To the authors' knowledge, no mapping-based approach has previously tackled the nonlinear spectra of such a large system.

II. BACKGROUND THEORY
We consider a molecular system consisting of a single electronic ground state and a manifold of coupled excited states. The Hamiltonian for such a system can be written as: 51,52Ĥ where both the nuclear configuration x and momentum p are vectors which have been mass weighted such that all degrees of freedom have unit mass. The ground-state potential operator isV g (x) = |0 V g (x) 0|. 53 The remaining potential operators describe all possible excited electronic states associated with a F chromophore system and are partitioned into their constituent subspaces according to the number of excitons present:V e (x) is the F ×F potential operator associated with the single-exciton subspace, |1 ,. . . ,|F , whileV ee (x) encodes the contributions from the double-exciton subspace, |1, 2 ,|1, 3 ,. . . ,|F − 1, F , of dimension 1 2 F (F − 1). We set = 1 throughout. All electronic couplings between the distinct exciton subspaces are neglected such that transitions are only induced when the system interacts with electromagnetic radiation. This field-matter interaction is accounted for by the dipole-moment operator,μ =μ + +μ − , which couples subspaces whose exciton numbers differ by one: whereâ † n (â n ) creates (annihilates) an exciton on chromophore n and µ n is the transition dipole moment associated with this chromophore.
Spectroscopic quantities of interest can then be calculated from Fourier transforms of the time-dependent polarization, which is induced in the medium by the applied radiation field. One common way to proceed is to perturbatively expand the time-dependent polarization in the field-matter interaction, which leads to the definition of the so-called optical response functions, where S (m) (t 1 , t 2 , · · · , t m ), is the the mth-order term in the expansion. 1 These optical response functions depend solely on the initial density matrix of the system, the dipole-moment operators [Eq. (2)] and the field-free dynamics associated with the Hamiltonian given in Eq. (1). We begin by considering the linear response in Sec. III before turning to the nonlinear case in Sec. IV.

III. LINEAR OPTICAL SPECTROSCOPY
The linear optical response can be described though the single-time correlation functions: where the key difference between these two functions is the initial state of the system. J abs (t) gives rise to the linear absorption spectrum, where the system is initialized in a thermal distribution associated with the electronic ground state, withρ g = e −βĤg andĤ g = |0 1 2 p 2 0| +V g (x). J fluor (t) gives rise to the fluorescence spectrum, where the system is initialized in the thermal state associated with the single-exciton subspace, witĥ ρ e = e −βĤe andĤ e = F λe=1 |λ e 1 2 p 2 λ e | +V e (x). Both of these correlation functions contain a forward (e −iĤet ) and a backward (e iĤgt ) propagator, which correspond to dynamics purely within the single-exciton subspace and electronic ground state respectively.
The linear absorption or fluorescence spectra are then obtained as follows: where θ(t) is the Heaviside step function and J(t) is the corresponding correlation function from one of Eqs. 3. The additional constants included in Eq. (4a) normalize the spectrum.
Quasi-classical expressions for single-time correlation functions such as Eqs. (3) can be easily obtained within both fully and partially linearized mapping-based techniques. These offer a practical way of calculating linear spectra for large condensed-phase systems based on phase-space averages over ensembles of classical trajectories. We will now briefly review these two approaches.

A. Fully Linearized Mapping Methods
Single-time correlation functions are expressed in a fully linearized mapping approach by representing the initial and final exciton state using the same single set of Cartesian mapping variables. We represent these as Z = {Z 0 , · · · , Z F }, a set of N = F + 1 complex numbers whose real and imaginary parts correspond to the usual MMST mapping variables. Thus the correlation functions given in Eqs. (3), are approximated by the phasespace averages: where µ + (Z) = 1 2 λe Z * λe µ λe Z 0 is the mapping representation of the dipole operatorμ + , |0 is the electronic ground state and |λ e is one of the single-exciton states. The dipole operator,μ − , at time t is then represented by the time-evolved mapping variables, µ − (Z(t)) = λeλ ′ e Z * λe λ e |ρ e (x, p)|λ ′ e µ λ ′ e Z 0 . The excitonic dynamics must be performed in the combined space of N = F + 1 components, Z, because the correlation functions are represented by a single set of mapping variables within the fully linearized approach. Therefore the dynamics associated with the forward and backward propagators, which occur in the single-exciton subspace and ground state respectively, cannot be treated separately. Additionally, the nuclear operators are described in terms of their Wigner transform. Note that, at each point in the nuclear phase-space,ρ e (x, p) is an F × F matrix in the single-exciton subspace, while ρ g (x, p) is a scalar function.
In Eq. (5), the phase-space average over the initial values for the Cartesian mapping variables, Z, and the nuclear coordinates, x and p, is given by: · · · = dx dp dZ · · · ρ m (Z), where ρ m (Z) is the mapping-variable distribution from which these variables are initially sampled. As outlined above, this fully linearized mapping-based approach for calculating optical absorption spectra has been also used in previous work. 11 In this paper, the fully linearized method we consider is spin-LSC, 36 where the excitonic degrees of freedom are represented using spin-mapping on the W-sphere. It has been already shown that such a method is able to reproduce single-time correlation functions to a high accuracy. 35 This is because spin-mapping, unlike for MMST mapping, exactly treats the identity operator within the underlying theory, leading to superior accuracy in calculating identity containing correlation functions (similar to the benefits seen in Refs. 28, 29, 31, and 54). Additionally, all configurations in the spin-mapping space correspond to a valid state in the real exciton space, alleviating the problems with MMST mapping associated with classical trajectories being able to 'leak' out of the physical mapping subspace. 55 For spin-LSC, the mapping-variable distribution corresponds to sampling the Cartesian mapping variables uniformly from the surface of a hypersphere of radius R 2 W = 2 √ N + 1, i.e., which is referred to as full-sphere sampling. The equations of motion for the mapping and nuclear phase-space variables then take the standard form, equivalent to those used in MMST mapping: whereV (x) =V g (x) +V e (x) is the potential matrix associated with the combined ground state and single-exciton subspace (i.e., a direct product) and is the mapping-variable representation of the force operator. Additionally, γ = (R 2 W − 2)/N is the generalized spin-mapping zero-point energy parameter, as derived in Ref. 36.
A key disadvantage of this fully linearized approach for calculating linear spectra is that for harmonic systems, the method is still not exact even when there are no diabatic couplings between the chromophores. However, in this case, even the much simpler WACL, 23 for which the nuclear force is associated with the time-independent arithmetic mean of the ground and a single excited state surface (i.e., F = − 1 2 (∇V g (x) + λ e |∇V e (x)|λ e ), where |λ e is one of the states in the single-exciton subspace), can exactly reproduce linear spectra. We would therefore like to design mapping-based methods which reduce to WACL in this special case. One way of fixing this problem is by 'quantizing' the mapping variables in terms of action-angle variables. This approach is employed by the optimized mean-trajectory (OMT) method, 41-43 which initializes half populations in both the ground and one of the excited states. An alternative is provided by the partially linearized approach, which can also be made to be exact in this limit as we describe in the following.

B. Partially Linearized Mapping Methods
The main difference between fully and partially linearized approaches is that in the latter, the electronic forward and backward propagators are represented using separate sets of mapping variables. The nuclear dynamics are however still described using a single set of nuclear variables. For the correlation functions given by Eqs.
(3) only the forward propagator corresponds to coupled exciton-nuclear dynamics within the single-exciton subspace, with the backward propagator involving only nuclear dynamics in the ground state. This means that only one set of mapping variables is actually needed to correctly describe the excitonic dynamics associated with these functions. Therefore expressions for the correlation functions within a partially linearized approach are given by: where Z e = {Z 1 , · · · , Z F } is now a set of Cartesian mapping variables purely within the single-exciton subspace with N = F components. The (partial) Wigner transform of the initial density matrix, ρ g (x, p) andρ e (x, p), and all phase-space averages over classical variables are defined equivalently as in the fully linearized approach, the key difference being that we use the smaller space, Z e . The time-evolved kernelŵ e (Z e , t) =Û e (t)ŵ e (Z e ) is a F × F matrix that is used to represent individual propagators associated with the single-exciton subspace and is evolved in time for each trajectory using the time-ordered propagator,Û e (t). At t = 0, the kernel matrix elements are defined by λ e |ŵ e (Z e )|λ ′ , where γ is the zero-point energy parameter associated with the mapping. In addition, the time-evolved propagator is defined as: U e (t) = e −iVeg(x(t))ǫ · · · e −iVeg(x(2ǫ))ǫ e −iVeg(x(ǫ))ǫ , (9) where ǫ is the time-step, and the potential matrix corresponds to the potential matrix associated with the single-exciton subspace defined relative to the ground-state potential. Defined in this way, the contribution of the time-evolved kernel to the partially linearized expressions for the correlation functions [Eq. (8)] also includes the e i t 0 dt1Vg(x(t1)) factor that arises from the quasi-classical expression for the backward ground-state propagator. For excited-state subspaces consisting of multiple excitons (required in Sec. IV), the associated potential matrices used for the exciton dynamics are also defined relative to the ground state potential such that this factor is removed when considering coherences which do not involve the ground state. The mapping variables are evolved under the standard equations of motion and the nuclear force is given as the average force associated with the forward and backward exciton paths: where V e (Z e , x) is the mapping-variable representation of the potential matrix associated with the single-exciton subspace. One of the features that we would like from a nonadiabatic dynamics approach for computing linear optical spectra is that it reduces to WACL in the absence of diabatic couplings between the chromophores and hence is exact for harmonic systems in this limit. The fact that the nuclear force used in WACL is also the average force associated with the single excited state and the ground state suggests that the partially linearized approach for calculating linear spectra can be devised in such a way that it fulfills this requirement. The key advantage of the partially linearized approach over WACL is that it includes the effects of nonadiabatic transitions and exciton energy transfer in systems when the chromophores are coupled.
In order to guarantee that the partially linearized approach does reduce to WACL in this case, the initial values for the mapping variables must be sampled using focused initial conditions in the single-exciton basis |λ e . This means that the mapping-variable distribution in Eq. (6) is given by: In this paper, we consider two partially linearized techniques both using focused initial conditions: the standard PLDM 45 approach, which corresponds to γ = 0 and the spin-PLDM 46,47 approach, which uses the same Wsphere spin-mapping approach as spin-LSC and therefore corresponds to γ = (R 2 W − 2)/N , where R 2 W = 2 √ N + 1 is the W-sphere radius. 36 Note, however, that because the size of the exciton space for computing linear spectra is different for spin-LSC and spin-PLDM (N = F + 1 and N = F respectively), the two methods actually use different values of the spin-radius and zero-point energy parameter. The advantage of treating the electronic ground state and single exciton subspace separately in spin-PLDM is that the most appropriate zero-point energy parameter for each distinct subspace can be used, which is expected to lead to more accurate results.
For standard PLDM, Eq. (8a) gives identical results to a previously implemented partially linearized approach for calculating optical absorption spectra. 49,50 However the implementation of our expression for the associated single-time correlation function has been simplified, so that only one set of mapping-variables is used for the simulation and the dimension of this set is reduced to the size of the single-exciton subspace (instead of considering the combined space of the electronic ground state and singleexciton subspace for both forward and backward paths as in the original approach). Formally, this has no effect on the accuracy of the standard PLDM method because the focused initial conditions treat the double combined space identically to the separated spaces. However, as we discussed above, it has important implications for spin-PLDM.
We showed in previous work that using focused initial conditions with spin-PLDM does not significantly affect the accuracy of obtained single-time correlation functions, but does mean that an order-of-magnitude fewer trajectories are needed to reach convergence of the mapping-variable integrals. 47 In contrast, using focused initial conditions with standard PLDM is known to rapidly degrade the quality of the results with increasing simulation time. 45 Therefore one advantage of using spin-PLDM over standard PLDM is that with focused initial conditions, the method gives accurate results for longtime dynamics when couplings exist between the chromophores, while still retaining the connection to WACL. In our previous work, it has also been been observed for a range of model systems that spin-PLDM generally exhibits greater accuracy when calculating single-time correlation functions compared to standard PLDM, even when focused conditions are not implemented. In particular, even though standard PLDM is reasonably accurate at short times, spin-PLDM is observed to be even better. 46 Because the linear optical response functions generally decay rapidly to zero, spin-PLDM is hence ideally suited to accurately calculate such quantities. We now apply these techniques to calculate linear spectra for excitonic condensed-phase systems.

C. Results
We consider exciton systems consisting of F chromophores, where each chromophore is linearly coupled to its own independent harmonic bath. The full Hamiltonian is given by: where f is the number of nuclear degrees of freedom associated with a single chromophore and ω shift shifts the chromophore energies such that F n=1 ǫ n = 0. 56 From this expression, the Hamiltonian associated with each of the various exciton subspaces can then be obtained. First, the ground-state Hamiltonian corresponds to the bath Hamiltonian projected into the zero exciton subspace:Ĥ g = |0 H B 0|. Additionally, the Hamiltonians associated with the other exciton subspaces then contain their corresponding projected bath Hamiltonian, along with an excitonic state matrix containing the associated matrix elements ofĤ S andĤ SB between all of the basis states of the subspace.
For all the exciton models considered in this paper, the distribution of the nuclear frequencies within each of the baths and their couplings is determined by the Debye spectral density: where Λ is the reorganization energy and ω c is the characteristic frequency of the bath. In order to have a finite number of nuclear degrees of freedom for our trajectory simulations, the continuous bath is discretized using the scheme employed in Ref. 57. When calculating linear spectra of nonadiabatic systems, quantum master equations are widely used. Of these, Redfield theory 52 is the most common, which is derived from a perturbative expansion in the excitonnuclear coupling and employs a Markovian approximation. In this paper, we always apply the secular approximation to Redfield theory, which has the advantage that the obtained dynamical populations are guaranteed to remain positive 52 and thus removes the possibility of obtaining unphysical negative peaks within calculated spectra. 20 A more sophisticated approach is the second-order time-convolutionless (TCL2) master equation, 14,58,59 which while still perturbative in nature is based on the second-order cumulant approximation such that it is by construction exact for systems with no diabatic couplings interacting with a harmonic bath. Such methods can be easily applied to excitonic systems and provide an interesting comparison with our partially linearized approach, which is the main focus of this paper.
In order to calculate the linear spectra using mappingbased methods, the (partial) Wigner transform of the Boltzmann operator for the appropriate excitonic subspace must first be calculated. For a harmonic bath, the Wigner transform of the Boltzmann operator in the exciton ground state is given by where α j = tanh( 1 2 βω j ). Due to the presence of excitonnuclear coupling, the Wigner transform of the Boltzmann operator in the single-exciton subspace is challenging to obtain exactly. 60 For simplicity, we approximate this quantity as follows: where this expression is obtained by taking the Wigner transform of the following approximate Trotter splitting for the single-exciton Boltzmann operator:ρ e ≈ e − β 2Ĥ S e −β(HB+ĤSB) e − β 2Ĥ S and is therefore accurate to second order in the diabatic couplings, ∆ nm and to first order in the exciton-nuclear coupling coefficients, c j . In principle, we could go beyond this approximation, but it is found to be sufficient for the cases we have tested. For the quantum master equations, the Boltzmann operator in the single-exciton subspace is approximated using only the purely excitonic Hamiltonian (i.e.,ρ e ≈ e −βĤS ) in accordance with previous work. 20 All the mapping-based methods considered here can exactly reproduce the time-dependent quantum dynamics of a bare electronic system. Also for a harmonic bath, our treatment of the bare nuclear dynamics by sampling the initial nuclear coordinates from the exact Wigner density and propagating them using classical equations of motion is also exact. Hence any errors appearing in the calculated optical spectra will solely arise due to the approximate description of the exciton-nuclear coupling. This therefore acts as a good test to compare the relative accuracies of calculating spectra with different mappingbased methods.

The Frenkel Biexciton Model
We first consider a previously studied Frenkel biexciton model, 11,12,14,42,44,49,50 which allows the accuracy of our partially linearized approach to be easily compared to other approaches. We use the same parameter sets considered in Ref. 11, which offer a comprehensive test over the various regimes that the model encompasses, from low-to high-temperature and homogeneous to inhomogeneous broadening. As in previous work, the transition dipole moments associated with the two exciton sites are chosen to be antiparallel with µ 1 = −5µ 2 . Additionally, exact results can be obtained using the HEOM approach, which we calculated using the open source pyrho package. 61 All results presented in this paper obtained with partially linearized approaches sample the mapping variables using focused initial conditions, for the reasons discussed above. This is illustrated in Fig. 1, where the optical absorption spectrum for two low-temperature (T = 72 K) Frenkel biexciton models with no diabatic coupling are reproduced exactly using partially linearized approaches either with the original MMST mapping (PLDM) or spinmapping (spin-PLDM). Because of the additional linearization approximation applied in the excitonic subspace, fully linearized approaches such as spin-LSC cannot be made exact in this case and exhibit large errors. Figure 1 also illustrates the advantage of using the TCL2 approach over Redfield, as only the former is exact in this case. From a theoretical perspective therefore, partially linearized approaches and TCL2 are better engineered to  correctly describe the dynamical coherences when calculating linear spectroscopic quantities of interest. Even though neither partially linearized approaches nor TCL2 can exactly reproduce the nonadiabatic dynamics of exciton models with diabatic couplings, the advantage of both these methods is also demonstrated by the reasonable accuracy of the numerical results when calculating optical absorption spectra. First we consider the relatively slow-bath and high-temperature models (i.e., small ω c and large T ), where the corresponding absorption spectra calculated using a range of methods are given by the two right-hand columns of Fig. 2. For both these models, all of the mapping-based methods can almost perfectly reproduce the absorption spectra, due to the following reasons. All such techniques are known to be able to reproduce the short-time dynamics essentially exactly in the high-temperature limit 46 and because all mapping-based methods exactly describe Rabi oscillations for an isolated excitonic system, they all by construction are able to correctly describe the static-nuclear limit, synonymous with a slow nuclear bath. The fact that a static-nuclear approximation (labelled 'Static' in Fig. 2) can also reproduce the spectrum extremely accurately, where the classical nuclear variables are still sampled from the Wigner distribution, ρ b (x, p), but are not evolved in time, confirms that the spectrum is almost entirely dominated by inhomogeneous broadening and hence these models do not correspond to a particularly challenging regime. Redfield theory, however, nonetheless completely fails in this case, as its Markovian approximation can only correctly describe the coupling to fast high-frequency modes. 11,14,59,62 In contrast to the relatively slow-bath models, the static-nuclear approximation cannot reproduce the absorption spectra for the relatively fast-bath and lowtemperature Frenkel biexciton models given by the two left-hand columns of Fig. 2, confirming that homogeneous broadening dominates here and that these models constitute a more challenging regime for accurately obtaining absorption spectra. Of the mapping-based approaches tested, spin-LSC is seen to produce the largest error,  in keeping with the fact that fully linearized approaches cannot correctly describe the dynamical coherences generated when measuring spectroscopic quantities and so cannot even reproduce the simple limit of the model without diabatic couplings, as discussed earlier in this section. Similar errors have also been observed for other fully linearized approaches when calculating absorption spectra for these same challenging fast-bath and low-temperature models. 11 TCL2 also outperforms Redfield for the same reason. Despite this, fully linearized methods are still expected to accurately reproduce linear spectra in the high-temperature and static-nuclear limits, where the associated dynamics are generally observed to be extremely accurate. 29,35,36,46 Although both the partially linearized approaches and TCL2 are not able to exactly reproduce the absorption spectra for these relatively fast-bath models given in Fig. 2, they are able to qualitatively reproduce the important features which is sufficient for most applications. Spin-PLDM is shown to be slightly more accurate than PLDM (although here the difference is not dramatic), which is in line with our previous observations that spin-PLDM gives consistently more accurate short-time dynamics compared to other mapping-based approaches. 46 Figure 3 gives the fluorescence spectra calculated for a range of different methods for the same Frenkel biexciton models. The quantum master equation approaches (i.e., Redfield and TCL2) are shown to produce much less accurate fluorescence spectra in comparison to the corresponding absorption spectra. This is because for quantum master equations the nuclear degrees of freedom have been integrated out, which means that the Wigner-transformed thermal Boltzmann operator associated with the single-exciton subspace,ρ e (x, p), cannot be accurately described and at best can only be approximated using the purely excitonic part of the Hamiltonian,Ĥ S . In contrast, as for the absorption spectra, all of the mapping-based approaches can almost exactly reproduce the fluorescence spectra for the relatively slow-bath and high-temperature models (two right-hand columns of Fig. 3), because both the dynamics and the approximate expression used for the Wigner-transformed thermal Boltzmann operator associated with the single-exciton subspace [Eq. (15)] are exact in the high-temperature and static-nuclear limit. The fact that the partially linearized approaches can also qualitatively reproduce the fluorescence spectra for the relatively fast-bath and low-temperature models (two left-hand columns of Fig. 3) suggests that Eq. (15) is also sufficient for calculating fluorescence spectra away from the high-temperature limit, at least in this case. Again, we observe that spin-PLDM is slightly more accurate at obtaining fluorescence spectra than PLDM.
Using the same methods and approach, the absorption, circular dichroism and fluorescence spectra have also been calculated for a seven-state FMO model 11,20 at a range of temperatures, the results of which are given in the supplementary material. Our partially linearized approach is again able to accurately reproduce the linear spectra associated with this system. While the differing accuracies of the various methods are not so pronounced when applied to the FMO model compared to the Frenkel biexciton models studied above, similar conclusions can nevertheless be drawn.

IV. NONLINEAR OPTICAL SPECTROSCOPY
In an analogous fashion to the linear case, a vast array of nonlinear spectra can be obtained by applying various Fourier transforms to multi-time response functions. 1 In many systems the second-order response vanishes, and we will thus focus on the third-order response function, S (3) (t 1 , t 2 , t 3 ). This is the basis for 2D optical spectroscopy, which is particularly popular for investigating exciton relaxation, dephasing and transfer processes. Nevertheless the theory presented here can in principle be applied to calculate any nonlinear spectrum of any order through its associated optical response function.
In order to calculate a 2D optical spectrum, the following contributions to the third-order optical response function are needed: 1,3,20,63 whereÂ(t) = e iĤtÂ e −iĤt is the time-dependent Heisenberg representation of an arbitrary operatorÂ. Each contribution can then be represented by a double-sided Feynman diagram (Fig. 4) which illustrates its associated dynamical pathway through the exciton subspaces. These diagrams can be generated from their associated mathematical expressions [Eqs. (16)] as follows. First, the propagators for each Heisenberg operator are written explicitly and any propagators directly next to each other are combined. The exciton subspace associated with each combined propagator is then determined by starting at the initial density matrix in the electronic ground statê ρ g and noting thatμ + (μ − ) adds (removes) an exciton.
From Fig. 4, we see that each of the contributions to the 2D optical spectrum has the following features in common. First, for the initial and final time intervals (t 1 and t 3 ) commonly referred to as the evolution and detection times, the system is in a coherence between two distinct exciton subspaces. Fourier transforms are therefore applied over these times to give the 2D spectrum. The so-called rephasing and nonrephasing diagrams (denoted 'RP' and 'NR') 3,63 correspond to the forward or backward propagator occupying the lower energy exciton subspace of the coherence during the evolution time, t 1 . Second, the middle time interval (t 2 ) commonly referred to as the delay time, is fixed for a given 2D optical spectrum. The associated dynamics occur within a single-exciton subspace, although the forward and backward paths may however occupy different exciton states. The difference between the so-called ground-state bleaching and stimulated emission diagrams (denoted 'GB' and 'SE' respectively) is that the ground state, |g , is occupied during the delay time for the former, while the single-exciton subspace, |e , is occupied for the later. In addition, the excited-state absorption diagrams (denoted 'ESA') differ from the stimulated emission diagrams by occupying a coherence state during the detection time which involves the double-exciton subspace, |ee . From these contributions, the 2D optical spectrum can be obtained by performing the following Fourier transform: 3,20 where Eqs. (17b) and (17c) are valid for both the rephasing and nonrephasing terms. Additionally, the impulsive pump-probe spectrum, also known as the transient absorption spectrum, corresponds to the case where the first two light-matter interactions coincide (i.e., when t 1 = 0) and hence the spectrum is given by: or equivalently using the nonrephasing diagrams.
In order to compute nonlinear spectra using mappingbased techniques, quasi-classical expressions for multitime correlation functions, such as those given by Eqs. (16), must be derived. For fully linearized methods it is not obvious how to accomplish this, because the phase-space integrals over the mapping-variable representations of dipole operators can only at most correctly reproduce the trace of the product of two operators. 35 The OMT method has however been used to compute multi-time correlation functions by incorporating discontinuous jumps in the mapping variables at the end of each time interval to mimic the transitions induced by the field-matter interaction. 41,42 While this offers a viable approach for computing nonlinear spectra using fullylinearized approaches, we note that it does not reduce to WACL when the system has no diabatic couplings between the chromophores. It was also found that for the pump-probe spectrum with t 2 = 0, the OMT approach can give different results for the GB and SE contributions to the full spectrum, even though the associated quantum correlation functions for each are identical in this limit. 42 Fully linearized mapping based techniques have additionally been used to obtain nonlinear spectra using a explicit light-field approach, where the excitonic mapping variables are evolved in time under the field-matter interaction. 12 This is in principle more general than the optical response function approach, and is thus able to compute spectra in strong electromagnetic fields. However, when the experiment is well described by the perturbation theory, it is cleaner to directly compute the nonlinear multi-time response functions. For example, using the explicit light-field approach to obtain the weakfield limit of the 2D optical spectra is more computationally expensive, because the simulation must be performed with twelve different pulse phases, in order to implement the phase-matching condition required in order to ob-tain the desired signal. Additionally, because the explicit light-field approach closely mirrors the experimental procedure, it is difficult to analyze the different contributions to the spectrum beyond what can be done from experiment. In contrast, because within the response function approach, each individual contribution to the 2D optical spectrum [Eqs. (16)] can be computed independently, the full signal can then be decomposed into the underlying parts associated with the corresponding processes.
The partially linearized approach overcomes many of the problems of the fully linearized approaches, in particular because it uses a different set of mapping variables to represent each propagator and can therefore be directly applied to calculating multi-time correlation functions. In Sec. IV A, we use such an approach to generate quasi-classical expressions for the multi-time correlation functions given in Eqs. (16).

A. Partially Linearized Mapping Methods
In order to generate partially linearized mapping based expressions for the third-order correlation functions, we employ an analogous approach to that outlined in Sec. III B for linear spectroscopy. We first consider the correlation functions which do not involve the doubleexciton subspace. As before, the propagators associated with the single-exciton subspace are each represented with a time-evolved kernel containing an independent set of mapping variables purely within this subspace. While the propagators associated with the ground state do not need to be represented by mapping variables, their effect is still accounted for by propagating the single-exciton subspace kernels using the associated potential matrix defined relative to the ground state potential, as described in Eq. (9). Finally, the nuclear force associated with each time interval is given as the average of the forces associated with the forward and backward exciton paths, as before. Hence using the Feynman diagrams associated with the contributions to the third-order optical response function (Fig. 4), which explicitly show the forward and backward exciton paths during each time interval, the partially linearized expressions can be easily generated. For example, the stimulated emission rephasing correlation function can be computed by: In practice, the expression for this three-time correlation function can be evaluated as follows. First, the nuclear phase-space variables are sampled from the initial Wigner density, ρ g (x, p) and the mapping variables Z e from focused initial conditions. All these classical variables, along with the associated Stratonovich-Weyl kernel,ŵ e (Z e ), are then propagated for the maximum considered t 1 time using the nuclear force given by the first line of Eq. (19b) and the values of these quantities are retained for each intermediate time-step. Second, a new set of mapping variables, Z ′ e , are sampled and a new kernel,ŵ e (Z ′ e ), is generated. For each intermediate t 1 time, the saved classical coordinates (x(t 1 ), p(t 1 ), Z e (t 1 ) and Z ′ e ) along with both kernels (ŵ e (Z e , t 1 ) and w e (Z ′ e )) are propagated for the t 2 delay time using the nuclear force given by the second line of Eq. (19b). Finally, the classical coordinates at the end of the t 2 timeinterval, along with the kernelŵ e (Z ′ e , t 2 ) are propagated for the maximum considered t 3 time using the nuclear force given by the final line of Eq. (19b) and the matrix elements of the kernel are retained at each intermediate time-step. The contribution to this correlation function [Eq. (19a)] is then calculated explicitly for each intermediate t 1 and t 3 using the corresponding time-evolved kernels. This data is then Fourier-transformed to generate the required spectrum.
For the excited-state absorption correlation functions, an additional time-evolved kernel containing mapping variables purely within the double-exciton subspace, w ee (Z ee ), is required. In the same way as before, the partially linearized expression for the excited-state absorption correlation functions and the associated nuclear forces can be generated from their associated Feynman diagrams (Fig. 4), which for the rephasing diagram gives rise to: The implementation for this correlation function is similar as for the stimulated emission correlation function considered previously, except that in general a third set of mapping variables Z ee is sampled after the t 2 delay time.
While a partially linearized approach has already been used to calculate third-order optical response functions, 50 our approach differs in the following ways. First, spinmapping instead of MMST mapping is used to describe the excitonic dynamics. Not only is spin-PLDM seen to offer an improvement over standard PLDM for linear spectra (Sec. III), it has also previously been observed to give rise to superior accuracy in obtaining population dynamics. 46 This difference is particularly pronounced when focused initial conditions are implemented, because they are known to significantly degrade the accuracy of standard PLDM. 45 On the other hand, using focused initial conditions in a partially linearized approach is particularly advantageous, as by construction the method will correctly reduce to WACL in the absence of exciton couplings within each distinct subspace. Focused spin-PLDM is therefore the method of choice, as it can simultaneously obtain correlation functions extremely accurately, while still retaining this connection to WACL. Second, we calculate each contribution to the third-order optical response function separately using its own specific mapping-based expression, which allows the dynamics in each of the distinct excitonic subspaces to be treated separately. While this makes no difference to the obtained results when using focused standard PLDM, the different dimensions of the considered exciton spaces for the two approaches will result in different values for the zeropoint energy parameters when using spin-PLDM. An advantage of our approach therefore is that by using the most appropriate zero-point energy for each of the distinct exciton subspaces with spin-PLDM should lead to more accurate results. Finally, in the original PLDM approach of Ref. 50, the mapping variables are resampled at the beginning of each time-interval, whereas in our approach the mapping variables are only resampled when a new exciton is created by a dipole operator. Because the resampling of mapping variables in a method generally makes the results more difficult to converge, we hence choose to do this only when it is absolutely necessary. Because of this difference, our standard PLDM approach for nonlinear spectra is in principle subtly different from that presented in Ref. 50, although in practice the obtained results are seen to be essentially identical for the systems treated in this paper.

B. Results
We now apply our partially linearized approach outlined in Sec. IV A to calculate pump-probe and 2D optical spectra. We study some of the same Frenkel biexciton models for which linear spectra were computed in Sec. III C as well as a seven-state FMO model. Specific details associated with these models, as well as implementation details of the various methods, can be found in Sec. III C and the supplementary material. Because all of these models are harmonic, exact results can again be obtained using the Hierarchical equations of motion (HEOM) approach, which we calculated using the opensource pyrho package. 61

The Frenkel Biexciton Model
An interesting aspect of the nonlinear spectra we consider is that they also probe the double-exciton subspace. For the Frenkel biexciton model, this subspace consists of a single state where both chromophores are occupied, such that mapping variables are not needed to describe the double-exciton dynamics. The corresponding Stratonovich-Weyl kernel in the partially linearized expressions for the ESA response functions [for example, Eq. (20a)] can thus be replaced as follows: which is essentially the single state analogue of Eq. (9). Figure 5 presents the pump-probe spectra for both a low-and high-temperature Frenkel biexciton model for various t 2 delay times. Only the methods from Sec. III C which are able to calculate nonlinear spectra within the optical response function approach are used. In particular, we do not present results of fully linearized methods as these cannot directly calculate multi-time correlation functions, although the reader may wish to compare our results with the fully linearized method of Refs. 11 and 12 obtained using the explicit light-field approach. For the relatively slow-bath and high-temperature model (i.e., small ω c and large T ), given by the bottom row of Fig. 5, the spin-PLDM results are essentially indistinguishable from the numerically exact results, in keeping with what was also observed for linear spectra. While the standard PLDM approach is also able to accurately reproduce the high-temperature pump-probe spectra at t 2 = 0, the error in the results increases for longer delay times. This is essentially consistent with earlier work, 46 where standard PLDM was found to exhibit significant errors when computing long-time population dynamics, even at high-temperature. In this case, even though standard PLDM still captures the correct qualitative features of the pump-probe spectra at non-zero delay times, the error is not much better than that exhibited by far simpler approaches, such as computing the exciton dynamics with static nuclei (labelled 'Static'). The fact that the static-nuclear approximation deviates from the exact results for these pump-probe spectra at non-zero t 2 times also illustrates that homogeneous broadening effects can still be present in nonlinear spectra even for models with relatively slow baths.
For pump-probe spectra associated with the relatively fast-bath and low-temperature Frenkel biexciton model, given by the top row of Fig. 5, all methods now exhibit errors compared to the numerically exact results. This is to be expected because for this relatively fast bath homogeneous broadening dominates, as illustrated by the fact that the static-nuclear approximation fails dramatically, and hence this system poses a greater challenge for accurately obtaining the associated nonlinear spectra. Spin-PLDM however is still able to qualitatively capture the correct features of the spectrum, even at large t 2 delay times and is significantly more accurate than standard PLDM.
Another advantage of the optical response function approach for calculating nonlinear spectra is that the full signal can be decomposed into its constituent parts associated with different dynamical pathways, giving greater Spin-PLDM PLDM Static Exact insight beyond what can be directly obtained from experiment or the explicit light-field approach. Figure 6 gives the three distinct contributions to the pump-probe spectra (GB, SE and ESA) for the two Frenkel-biexciton models considered in Fig. 5 at various t 2 delay times, calculated using spin-PLDM (coloured lines) and the numerically exact HEOM approach (dashed black lines). As discussed before, the spin-PLDM approach is able to qualitatively reproduce the main features of the various contributions to and thus the full pump-probe spectra in both the high-temperature and low-temperature regimes.
Considering each of these contributions in turn allows the full pump-probe spectra given in Fig. 5 to be better understood. First, the GB signal corresponds to the linear absorption spectrum for all t 2 delay times, because the initial nuclear density matrix,ρ g , is conserved by the ground state dynamics during the t 2 delay time. Second, the SE signal becomes the fluorescence spectrum in the limit of t 2 → ∞, as the system relaxes toρ e due to the dynamics in the single-exciton subspace during the t 2 delay time. This explains why the lowest-energy peak within the SE signal increases in intensity with increasing t 2 delay time, in accordance with Kasha's rule. Finally, the highest-energy peak in the ESA signal increases in intensity with increasing t 2 delay time, because absorption to the double-exciton subspace occurs from the remaining unoccupied state within the single-exciton subspace.
The 2D optical spectrum can also be calculated using the same approach. This is more computationally demanding than pump-probe spectra, as now the contributions to the optical response function must be computed on a (t 1 ,t 3 ) two-dimensional grid. The 2D optical spectra for the relatively slow-bath and high-temperature Frenkel biexciton model for various t 2 delay times is given in the supplementary material, which is the main benchmark used in previous work for nonlinear optical spectra. 12,14,42,50 As for the pump-probe spectra, the spin-PLDM 2D optical spectra for this regime are essentially indistinguishable from the exact results, while even though the standard PLDM results are qualitatively accurate, they do exhibit small but noticeable errors for large t 2 delay times. We however note that although the static-nuclear approximation does deviate from the exact results for this parameter set, it is nevertheless able to qualitatively reproduce the important features of the spectra, illustrating that this model does not pose a challenging test for our method. We therefore consider the 2D optical spectra for a more challenging relatively fastbath and low-temperature regime of the model, where homogeneous broadening effects dominate and the staticnuclear approximation fails. In this regime, shown in Fig. 7, both partially linearized approaches are able to qualitatively reproduce the most important features. However, for non-zero t 2 delay times, spin-PLDM is again more accurate than standard PLDM, in particular for the t 2 = 200 fs case.

The Fenna-Matthews-Olsen complex
An example of a biologically relevant system for which it is interesting to compute nonlinear spectra is the sevenstate FMO model. All the Hamiltonian parameters for this model can be found in Ref. 20. This system contains F = 7 single-exciton states and F (F − 1)/2 = 21 doubleexciton states and hence requires a set of mapping variables to describe the associated exciton dynamics within each of these subspaces. The excitonic space of this system is thus much larger than that of the Frenkel biexciton model considered previously (which contains two singleexciton states and only one double-exciton state), therefore acting as a good test for the viability of these methods to calculate nonlinear spectra in realistic condensedphase systems. To compute the pump-probe spectra, we assign dipole-moment orientations to each chromophore following Ref. 20 and consider the case where all laser pulses are equally polarized, such that the rotational averaging of these spectra can be performed by averaging the multi-time correlation functions over 10 representative electric-field directions (corresponding to the vertices of a dodecahedron) as described in Ref. 18. This averaging can be performed as a post-processing step after each trajectory is computed, to determine its contribution to the ensemble average. The same approach could also be applied to more complicated polarization sequences, where at most 21 electric field directions have to be considered. 64 Figure 8 gives the pump-probe spectra for the FMO model at T = 300 K, for three different t 2 delay times. As for the other high-temperature models, the spin-PLDM results are essentially indistinguishable from the benchmark, for both the components and the full spectra, with any errors being numerical in nature rather than systematic. As was also found for the Frenkel biexciton models, standard PLDM is accurate for the high-temperature pump-probe spectra at t 2 = 0 fs, but the error associated with the method increases for longer t 2 delay times. Spin-PLDM therefore clearly exhibits superior accuracy in computing nonlinear optical spectra than standard PLDM.
It seems that the only disadvantage of the spin-PLDM approach relative to standard PLDM is that for large exciton subspaces, the method requires a greater number of trajectories to converge. For the relatively small exciton subspaces of the Frenkel-biexciton models, there is no issue and both standard PLDM and spin-PLDM require about 10 6 trajectories to reach full convergence on the scale of the plot when calculating the multi-time correlation functions needed for nonlinear spectra. However when calculating the excited-state absorption correlation functions of the FMO model, for which the very large 21 state double-exciton subspace must be considered, we used 1 × 10 9 spin-PLDM trajectories to reach convergence compared to 2 × 10 7 trajectories when using standard PLDM. This is in contrast to the 1 × 10 7 and 4 × 10 7 spin-PLDM trajectories used to obtain the ground-state bleaching and stimulated emission correlation functions respectively and the 1 × 10 6 trajectories needed to obtain the linear absorption and fluorescence spectra.
We have, however, found a simple way to dramatically reduce the cost of the spin-PLDM approach for cal- culating these excited-state absorption correlation functions, while still retaining the superior accuracy of the method. For the double-exciton subspace, we choose to use γ = 0 for the focused initial conditions of the mapping variables, Z ee , the associated kernel,ŵ ee (Z ee ) and the mapping variable expression for the force operator, ∇V ee (Z ee , x), while the spin-PLDM expression for the zero-point energy parameter is still used for the singleexciton subspace. This essentially corresponds to a standard PLDM treatment of the double-exciton subspace, which is a reasonable approximation because the correlation functions decay rapidly as a function of the t 3 time and it is known that standard PLDM is also accurate for relatively short propagation times. For example, a maximum t 3 time of only 500 fs was required to correctly obtain the pump-probe spectra for this FMO model at T = 300 K. Hence treating the largest doubleexciton subspace using standard PLDM should drastically reduce the number of trajectories needed for convergence, without significantly affecting the accuracy of the method. Fig. S4 in the supplementary material gives the pump-probe spectra for our modified spin-PLDM approach, where we only used 4 × 10 7 trajectories to obtain the excited-state absorption correlation function. The accuracy of these results can be seen to be almost equal to the full spin-PLDM method, but now obtained at a significantly reduced computational cost similar to the standard PLDM approach. This significant improvement in efficiently should make it possible to compute the 2D spectra of FMO with our partially linearized approach, which requires additional scanning over the t 1 time. We will address this in future work.

V. CONCLUSIONS
In this paper we have shown how classical trajectory mapping-based methods can be used to compute optical spectra for nonadiabatic systems through the constituent response functions. In particular we use a partially linearized approach, which represents each individual propagator within any correlation function using independent sets of mapping variables. This means that single-and multi-time correlation functions, used to compute linear and nonlinear spectra respectively, can all be calculated on an equal footing. We have demonstrated that this approach can accurately and consistently reproduce the important features of absorption, fluorescence, pump-probe and 2D optical spectra for a range of systems and parameter regimes.
One way our approach differs from others is by requiring that the expressions for the correlation functions reduce to WACL when there are no diabatic couplings between the chromophores, such that the method is exact in this limit for systems containing harmonic nuclear degrees of freedom. While such a requirement can be satisfied using a partially linearized approach with focused initial conditions, it cannot be satisfied with fully linearized methods, often resulting in inaccurate linear spectra for systems in more extreme and challenging parameter regimes and generally giving rise to overbroadened peaks. Another distinction of our approach is that we treat the excitonic dynamics separately within the distinct subspaces (something which also cannot be done within a fully linearized approach); this eliminates unphysical forces which arise from mapping-variable contributions to exciton configurations that cannot occur and also means that the most appropriate spin-mapping zeropoint energy parameter for each subspace can be used, which is expected to result in more accurate dynamics. Overall, these points illustrate how partially linearized approaches are much more suited to computing spectroscopic quantities than their fully linearized counterparts.
We have considered two partially linearized mappingbased methods, standard PLDM and spin-PLDM, and have demonstrated that spin-PLDM consistently gives rise to the most accurate spectra. For high-temperature systems, the spin-PLDM results are essentially numerically exact, while we found that the method can still qualitatively reproduce the important spectroscopic features in more challenging low-temperature regimes. In contrast, standard PLDM was shown to exhibit errors in the nonlinear optical spectra at non-zero delay times, even in the high-temperature regime. One of the reasons that standard PLDM is less accurate than spin-PLDM is because focused initial conditions in standard PLDM are known to rapidly degrade the quality of the results, whereas focused spin-PLDM is observed to be just as accurate as the full-sphere sampling variant. 35,36,47 Therefore one of the advantages of spin-PLDM over standard PLDM is that focused initial conditions can be used to ensure that the method reduces to WACL in the absence of diabatic couplings, without compromising the accuracy of the obtained results for systems with coupled chromophores.
The only drawback of the spin-mapping version of PLDM is that the method is computationally expensive for systems with large exciton subspaces. To alleviate this problem, we introduced a modified spin-PLDM approach which treats the double-exciton subspace using the standard PLDM approach. Because the double-exciton subspace dynamics only occurs during the t 3 time, over which the multi-time correlation function rapidly decays, the relatively good short-time accuracy of standard PLDM should be sufficient to still accurately describe this. This modification was found to significantly reduce the number of trajectories needed to reach convergence, such that the new approach has a comparable computational cost to standard PLDM, while still retaining the superior accuracy of the original spin-PLDM approach. This therefore offers a practical way of computing nonlinear spectra in light-harvesting systems.
All our approaches are built on the perturbative response of the exciton-nuclear dynamics due to a weak coupling to the light field. However, for experiments involving intense laser fields, the effect of the classical electromagnetic radiation on the coupled exciton-nuclear dynamics must be included explicitly via a time-dependent Hamiltonian. In principle, it it will be possible to develop a partially linearized approach for this scenario which reduces to that already outlined in this paper in the weak-field limit. Because of this requirement, such a method will therefore differ from other previously considered mapping-based approaches to this problem. [11][12][13] Additionally for experiments involving extremely low intensity radiation, such as those performed within a cavity, a quantum description of light is needed to correctly describe the induced exciton-nuclear dynamics. This could be achieved within a partially linearized approach by treating each cavity mode as an additional harmonic degree of freedom within the system and then treating the dynamics of these modes classically. 65 Because spin-PLDM can in principle be applied to anharmonic problems, it also offers a route to computing linear and nonlinear spectra for more complex systems which cannot be treated by HEOM. Hence the importance of anharmonic effects on excitonic energy-transfer in photosynthetic light-harvesting systems could be investigated by performing atomistic simulations based on classical force fields, 66 leading to a greater overall understanding of the mechanism that gives rise to such highly efficient energy transfer in light-harvesting systems.

SUPPLEMENTARY
See the supplementary material for further results and for the data files of calculated optical response functions for spin-LSC, standard PLDM and spin-PLDM, which can be used to generate all the associated spectra presented in this paper.

ACKNOWLEDGMENTS
The authors would like to acknowledge the support from the Swiss National Science Foundation through the NCCR MUST (Molecular Ultrafast Science and Technology) Network. We also thank Graziano Amati, Johan Runeson and Joseph Lawrence for helpful discussions and comments.

DATA AVAILABILITY
The data that supports the findings of this study are available within the article and its supplementary material.
In this supplementary material, we present additional results for both linear and nonlinear spectra, which were not included within the main paper. While these results give rise to the same overall conclusions with regards to the relative accuracy of each of the considered methods, the discrepancies between the obtained results are not as large. They are nevertheless included here for completeness. We also present the results for a modified spin-PLDM approach for calculating the excitedstate absorption correlation functions used to compute 2D and pump-probe spectra. Through application to the Fenna-Matthews-Olsen (FMO) complex, we show that this new approach is just as accurate as the full spin-PLDM method, but requires significantly fewer trajectories to reach convergence. Data files containing all of the response functions corresponding to the spectra presented in this paper for standard PLDM, spin-LSC and spin-PLDM are also included. We also provide python scripts which can be used to perform the necessary Fourier transforms to generate the associated spectra.

I. CIRCULAR DICHROISM SPECTRA
In addition to electronic absorption and fluorescence spectra, an array of other linear spectroscopic quantities exist, which can be obtained as the Fourier transform of a single-time optical response function. An example is the circular dichroism spectrum, which measures the difference in absorption of right-and left-circularly polarized light. The associated correlation function is defined as: 1,2 J CD (t) = Tr e iĤgtm− e −iĤetμ+ρ g , (S1a) wherem =m + +m − andμ =μ + +μ − are the magnetic and electric dipole operators respectively andâ n is the exciton destruction operator for chromophore n. Addi-tionallyĤ g andĤ e are the Hamiltonians associated with a) Electronic mail: jonathan.mannouch@phys.chem.ethz.ch b) Electronic mail: jeremy.richardson@phys.chem.ethz.ch the ground and single-exciton subspaces, R n is the position vector associated with the centre of chromophore n and × corresponds to the vector cross-product. From these expressions, the fully linearized and partially linearized mapping expressions for the circular dichroism optical response function can be obtained in accordance with the expressions outlined in the main paper for absorption spectra. Additionally, the rotational averaging of the linear spectra over all possible electric field orientations can be achieved by considering the components of the two observable operators along each of the three Cartesian axes. 3 This is implemented as a postprocessing step after each trajectory is computed to determine its contribution to the ensemble average. In order to calculate the circular dichroism spectrum, the Fourier transform of the optical response function is taken as follows: where the positive part of the spectrum is normalized for positive frequencies. We now analyse the electronic absorption, circular dichroism and fluorescence spectra for the FMO complex at a range of temperatures, calculated using the quantum master equation and mapping-based methods considered in the main paper.

II. THE FENNA-MATTHEWS-OLSEN COMPLEX
For the FMO complex, we consider the same model as in Ref. 2, where complete details of the model and all the required parameters can be found. For all of the classical trajectory techniques, f = 60 nuclear degrees of freedom per chromophore were used. The numerically exact results were computed using the hierarchical equations of motion (HEOM); for the T = 30 K linear spectra, these were obtained from Ref. 2 while the other spectra were calculated using the open-source pyrho package. 4

A. Linear Spectra
The associated electronic absorption, circular dichroism and fluorescence spectra are given in Figs. S1, S2 and S3 respectively. Considering that the error of the static nuclear approximation is large in comparison with  the numerically exact results suggests that homogeneous broadening dominates. Such a regime is thus challenging and requires a dynamical theory for the nuclei in order to correctly reproduce the qualitative features of the spectra.
For the mapping-based methods (given in the first row of all three figures), we find that all of the generated spectra in the high-temperature limit (T = 300 K) are essentially indistinguishable from the exact results. This is in agreement with previous work, which has found that for the standard partially linearized density matrix (PLDM), spin-PLDM and spin-LSC methods, the relatively short-time dynamics become exact for hightemperature harmonic models. [6][7][8] However for the circular dichroism spectrum (Fig. S2) at this temperature, the standard PLDM result does exhibit a noticeably greater error than spin-PLDM.
As the temperature is decreased, the peaks associated with the spin-LSC linear spectra become more and more over-broadened compared to all of the numerically exact results. This is in agreement with the findings in the main paper, where this phenomenon was attributed to the fact that the fully linearized methods are unable to correctly describe the dynamical coherence between the ground and single-exciton subspaces. Although also not exact, the partially linearized results (standard and spin-PLDM) achieve closer agreement to the numerically exact spectra at low temperatures. In general though for all linear spectra and parameter regimes considered, spin-PLDM is found to be just as good or even better than standard PLDM at correctly reproducing the important spectral features.
For the quantum master equations (given in the second row of all three figures), the second-order timeconvolutionless (TCL2) quantum master equation is by far the most accurate and is able to almost exactly reproduce the electronic absorption and circular dichroism spectra at all temperatures. This was explained in the main paper by noting that the TCL2 approach is also able to correctly describe the dynamical coherence between the ground and single-exciton subspaces, as demonstrated by the fact that as for the partially linearized approaches, the method is exact for harmonic models without diabatic couplings. 9 This is not true for Redfield theory, which although seemingly accurate at low temperatures, exhibits large errors in the hightemperature limit, as has been noted in previous work. 2 In general, all quantum master equations are unable to correctly reproduce electronic fluorescence spectra because the initial coupled exciton-nuclear state cannot be correctly described by reduced densities for which the nuclear degrees of freedom have been integrated out.

B. Pump-Probe spectra
Here we present the pump-probe spectra associated with the FMO complex at T = 300 K, calculated using our modified spin-PLDM approach described in the main paper. Within this approach the ground-state bleaching (GB) and stimulated emission (SE) contributions to the spectra are calculated using the original spin-PLDM approach and are hence identical to the results presented in the main paper. The calculation of the excited state absorption (ESA) contribution to the spectra is however different, because the mapping variables in the doubleexciton subspace (Z ee ) are now treated with standard PLDM (γ = 0), while the mapping variables in the singleexciton subspace (Z e and Z ′ e ) are treated with the same zero-point energy parameter used in the original spin-PLDM method. This modification requires significantly fewer trajectories to reach convergence, because the standard PLDM focused initial conditions are more efficient than those for spin-PLDM, especially when sampling the large double-exciton space. The reason for this is explained in Ref. 12. By comparing the modified (Fig. S4), original spin-PLDM (given by Fig. 8 in the main paper) and standard PLDM results (Fig. S4), we see that the superior accuracy of the spin-PLDM approach is still retained in the modified approach, as desired.

III. THE FRENKEL BIEXCITON MODEL
The Frenkel biexciton model is a two chromophore system commonly used to benchmark methods used to compute electronic spectroscopic quantities. 5,[9][10][11]13 In particular, the relatively slow-bath and high-temperature version of the model (T = 300 K) has been used numerous times as a test for the calculation of 2D electronic spectra. [9][10][11] Even though this regime is less challenging than the relatively fast-bath and low-temperature case considered in the main paper, we include this here for completeness and in order to aid comparison with other results published in the literature. Further details of the model, including all the required parameters, can be found in Ref. 5. For all of the classical trajectory techniques, f = 100 nuclear degrees of freedom per chromophore were used. Exact results can also be computed using HEOM, where we have used the open source pyrho package. 4 The spin-PLDM 2D electronic spectra (given by the second row in Fig. S5) are essentially indistinguishable from the exact results (given by the first row in Fig. S5). This is again to be expected, as the dynamics produced by spin-PLDM have been observed to be numerically exact for high-temperature harmonic models. 8 While the standard PLDM results (given by the third row in Fig. S5) are also practically perfect for t 2 = 0, a small error can clearly be seen when t 2 is increased. This is most noticeable for the t 2 = 600 fs results, where even though standard PLDM is still able to qualitatively reproduce the correct features of the spectrum, the relative peak heights are different in comparison to the exact results. This error arises because standard PLDM is known to often be unable to correctly reproduce the long-time behaviour of dynamical observables, even in the hightemperature limit. 8