On the Applicability of the Caldeira-Leggett Model to Condensed Phase Vibrational Spectroscopy
Formulating a rigorous system-bath partitioning approach remains an open issue. In this context the famous Caldeira-Leggett model that enables quantum and classical treatment of Brownian motion on equal footing has enjoyed popularity. Although this model is by any means a useful theoretical tool, its ability to describe anharmonic dynamics of real systems is often taken for granted. In this Letter we show that the mapping between a molecular system under study and the model cannot be established in a self-consistent way, unless the system part of the potential is taken effectively harmonic. Mathematically, this implies that the mapping is not invertible. This ‘invertibility problem’ is not dependent on the peculiarities of particular molecular systems under study and is rooted in the anharmonicity of the system part of the potential.
pacs:78.20.Bh, 78.30.-j, 33.20.Ea, 05.40.Jc
Introduction and Theory.
Model systems play an important role for our understanding of complex many-body dynamics. Reducing the description to a few parameters can not only ease the interpretation, but enable the identification of key properties may11 (). In condensed phase dynamics the spin-boson leggett87_1 () and Caldeira-Leggett (CL) Caldeira-PRL-1981 (); Caldeira-AP-1983 () models have been the conceptual backbone of countless studies grabert88_115 (). The latter has been extended to describe linear and non-linear spectroscopy within the second-order cumulant approximation, termed multi-mode Brownian oscillator model in this context Mukamel1995 (); Tanimura-PRE-1993 (). It has become a popular tool for assigning spectroscopic signals in the last two decades.
The CL model comprises an arbitrary system (coordinates , potential ), which is bi-linearly coupled to a bath of harmonic oscillators (coordinates , potential ) via a system-bath coupling potential Caldeira-PRL-1981 (). Later Caldeira and Leggett extended the model to an arbitrary function of system coordinates in the coupling and motivated the linearity of the coupling on the bath side Caldeira-AP-1983 (). Here, we limit ourselves to the bi-linear version for the reasons that will become apparent later. Restricting ourselves to a one-dimensional case yields grabert88_115 ()
with the bath frequencies , unit bath masses, and the coupling strengths . Often the square on the right hand side of Eq. (1) is completed causing thereby a harmonic correction to the system potential : . We stress, that omitting this harmonic correction is not self-consistent, see Supplement. The CL model imposes no restrictions on the form of the system potential, apart from the presence of the aforementioned trivial correction.
The advantage of this model is the possibility to derive from it a simple reduced equation of motion termed generalized Langevin equation (GLE), that reads
where the bath is reduced to non-Markovian dissipation and fluctuations represented by the memory kernel and the stochastic force , respectively. The former is connected via a cosine Fourier transform to the spectral density
and the latter is a Gaussian process with zero-mean. Importantly, the derivation can be performed both in classical Mukamel1995 (); ZwanzigBook2001 () and quantum Ford-JSP-1987 (); Ford-PRA-1988 () domains. For the present purpose we would limit ourselves to the classical GLE, whose derivation from the CL model requires no further approximations. In this case, the fluctuating force is connected to the memory kernel by means of the so-called fluctuation-dissipation theorem (FDT)
which establishes the canonical ensemble with the temperature . We note that according to Ford et al. Ford-JSP-1987 (): “…if there is a universal description, then it must be of the form we have obtained [that is Eq. (2)]”.
The tremendous simplicity of the CL model is, thus, that it reduces the system-bath interaction, to a single spectral density. However, a real system Hamiltonian might very well differ from the CL one. In practical applications one, therefore, tries to find a mapping between the two in order to make full use of the model’s simplicity. Importantly, if such a mapping can be established, then the full quantum-mechanical treatment of the bath can be performed analytically via the Feynman-Vernon influence functional Feynman-AP-1963 (); Caldeira-PhysicaA-1983 () without any further approximations. Moreover, numerically exact hierarchy type equation of motion approaches are essentially based on the CL model Dijkstra-PRL-2010 (); Suess:2014gz (). Further, it is ideally suited for numerical methods that solve the Schrödinger equation in many dimensions giese06_211 (). The CL model can also be taken as a starting point for quantum-classical hybrid simulations, that is by treating only the usually low-dimensional system part quantum-mechanically and the bath (semi)classically egorov99_5238 (). Finally, the machinery for a purely classical treatment by means of a GLE is provided by the method of colored noise thermostats Ceriotti2009a (); Ceriotti2010 (). It is thus the CL model that establishes a unified framework for a quantum-classical comparison of dynamical properties in condensed phase systems provided that the real system can be mapped onto this model. We note in passing that various models based on a coupling that is non-linear in system coordinates exist, such as, e.g., the so-called square-linear CL model kuhn03:2155 (); Tanimura2009 (). They, however, correspond to different GLE types, which do not lead to any known mapping, see Supplement, and are, thus, not considered here. This puts forward the central question of this Letter, i.e. whether such a mapping onto the bi-linear CL model exists for an arbitrary system. Surprisingly, the amount of information concerning the applicability of this tremendously popular model is scarce. Makri argued for the applicability if the system-bath coupling can be treated in the linear response regime makri99:2823 (). Goodyear and Stratt derived the model in the short time limit introducing so-called instantaneous normal modes Goodyear1996 (). However, the bath statistics turned out to deviate from the predictions of the CL model at longer times Goodyear1996 (); Tuckerman1993a ().
In order to approach the central question, we discuss three possible ways for constructing a mapping onto the CL model. First, one might assume that the real bath can be either exactly or approximately described by the CL Hamiltonian in a certain coordinate system. This case corresponds to a direct representation of the bath in a harmonic form, Eq. (1). The system potential would be then the actual system potential as given, for instance, by an explicit force field; note the presence of the aforementioned harmonic correction in . Another possibility is to rigorously derive the GLE from first principles for an arbitrary open system. If the resulting equation coincided with that derived from the CL model, it would then serve as an indirect proof of the model itself. Here, Mori’s projection operator techniques are commonly involved to derive a GLE-type equation from the general Hamilton equations of motion, expressing noise and dissipation formally as projected quantities Mori1965 (); ZwanzigBook2001 (); Cubero-JCP-2005 (). The projection can be performed either in a linear or a non-linear way. The second mapping is based on the linear projection, which yields the GLE Eq. (2) and the FDT Eq. (4) without any further approximations. The system potential becomes effectively harmonic, i.e. at the price of projecting all the system anharmonicity into the bath. The effective harmonic frequency follows as with denoting canonical averaging. The third route employs a non-linear projection. Here, the system potential is given by the potential of the mean force, , felt by the system at coordinate , averaged over all bath coordinates Kawai2011 (); ZwanzigBook2001 (); Lange2006 (). In this formalism the memory kernel carries a functional dependence on the system’s positions and momenta and the fluctuating force has more complicated statistical properties than given by the FDT in Eq. (4) Kawai2011 (). This functional dependence is hardly tractable and an apparent simplification is to linearize the kernel with respect to positions and momenta. In this case the contribution proportional to vanishes and the resulting integral coincides with that in Eq. (2) Kawai2011 ().
These three routes provide us with the three versions of the GLE that are parametrized from explicit MD simulations upon inverting the equation
to obtain the memory kernel , where and are the momentum–momentum and momentum–system force time-correlation functions, see Supplement. Thus, the mapping between the explicit system under study and the CL model can be formulated as the mapping . The particular form of the system force, that is, the particular form of , specifies the GLE employed and thus the mapping , referred to as , and for the direct mapping, mapping based on linear and on non-linear projection technique, respectively. We note again that the latter mapping requires a linearization of the memory kernel. In order to verify a mapping and thus to deduce the applicability of the CL model, we compute the linear vibrational spectra from the GLEs and compare them against their explicit MD counterparts. The match between the spectra describes thus the success of the system-bath model, see Fig. S5 in the Supplement.
Models and Methods.
We have employed three hydrogen-bonded systems: HOD in HO, OH in HO, and the ionic liquid at room temperature. These examples have been chosen since H-bonding is known to come along with a considerable anharmonicity in the potential energy surface (PES). In particular, HOD in HO or DO is among the most studied hydrogen-bonded systems stenger01_027401 () and the comparison with the OH ion is used to elucidate the importance of the intramolecular forces. The ionic liquid represents a system, where, in addition to moderate H-bonding, there exists a strong Coulomb interaction between the ion pairs roth12_105026 ().
The HOD/OH in HO systems are comprised of 466 molecules in a periodic box with the length of nm interacting according to the force field taken from Ref. Paesani-JCP-2010 (). In the latter case the HOD molecule is substituted by the diatomic OH fragment. The ionic liquid simulations were carried out in a periodic box with the length of nm containing 216 ionic pairs and the force field described in Ref. Koeddermann2007 (). All explicit molecular dynamics (MD) simulations were performed with the GROMACS program package (Version 4.6.5) GROMACS (). The spectra were computed for the OH stretch in the HOD/OH molecules and for the CH stretch of the imidazolium ring Koeddermann2007 () in the ionic liquid, see Supplement for details. In the latter case the potential for the C-H stretching motion was re-parametrized to a Morse potential using DFT-B3LYP calculations roth12_105026 (). Note that in the spirit of the system-bath treatment the O-H bond lengths were taken as the respective system coordinates, , whereas all other degrees of freedom constituted the bath. We employed the “standard protocol” for calculating IR spectra, that is, a set of NVE trajectories, each ps long (time step fs), was started from uncorrelated initial conditions sampled from an NVT ensemble. The functions for the system coordinate were Fourier-transformed to yield the vibrational spectra Ivanov-PCCP-2013 (); note that these spectra are proportional to the absorption spectra in the one-dimensional case. In order to achieve convergence, 1000 trajectories for the considered stretching motion were employed, see Supplement for further details. As a numerical solver for the GLE we adopted the method of colored noise thermostats Ceriotti2009a (); Ceriotti2010 (). The simulation setup for GLE simulations was the same as that for the explicit ones.
Following the three routes described above, we start with checking whether the direct mapping exists, which implies that the actual Hamiltonian of the system studied has the form of Eq. (1). The linearity of the coupling on the system side can be verified via explicit MD simulations in a straightforward manner. Here, we sampled uncorrelated configurations and varied the system coordinate, , within the range accessible due to its thermal fluctuations keeping all bath coordinates fixed, see Supplement. The system-bath coupling probed this way was least-squares fitted to linear functions . In order to quantify the fit error with respect to a reasonable scale, we considered the relative deviation
with being the initial/final values of the probed range, respectively. The results averaged over the chosen dataset are depicted in the upper panels of Fig. 1 with the error bars given in red. It becomes obvious that excellent fit results can be obtained for the two aqueous systems (panels a and b) where the relative deviations are below %. In contrast, the averaged fit errors for the ionic liquid (panel c) are in the range of on the scale of the overall change of the potential on the accessible interval. Hence, we can conclude that for aqueous systems the coupling is linearly dependent on the bond length, , as it is assumed in the CL model, whereas it is by no means true for the ionic liquid considered.
Given the fact that the CL model is strictly equivalent to the GLE, the assumptions concerning the nonlinearity of the coupling on the bath side and the harmonicity of the bath, being hardly separable, together imply Gaussian statistics of the noise. The Gaussian statistics can be easily probed via explicit MD simulations by calculating the distribution of the environmental forces exerted on the system coordinate, which is kept fixed to remove the dissipative contributions, see Supplement. The resulting distributions shown in the lower panels of Fig. 1 are clearly asymmetric and centered at positive values for the two aqueous systems (panels d and e). The distribution for the ionic liquid (panel f) is centered almost at zero and is only slightly skewed, thus looking much more similar to a Gaussian. To summarize, checking the linearity on the system side suggests that the model should perform well for aqueous systems and badly for the ionic liquid, whereas the corresponding results for the linearity on the bath side imply the opposite. Another test for the bath’s linearity is provided by probing the temperature dependence of the spectral density, which should be absent if the two are strictly fulfilled. The results lead to similar conclusions, see Fig. S7 in the Supplement.
In order to elucidate the impact of these non-linearities onto observables we compare the GLE spectra with the ones obtained from explicit MD simulations in Fig. 2. The GLE spectra (green curves therein) for aqueous systems (panels a,b) are significantly different from the reference ones (black curves therein) both in shape and position. For the ionic liquid (panel c) these differences are not so pronounced, though still noticeable. This comparison supports the conjecture that the Gaussian statistics of the bath is more important, although it cannot be considered as its proof in general. Still, in all cases considered, the match of spectra is not achieved and thus the direct mapping fails to exist for the systems studied.
The second route is based on the mapping , which is exact and therefore the corresponding spectra (red stars in Fig. 2) perfectly match the explicit MD ones (black lines) for all systems considered. Importantly, this mapping is established at the price of putting all the system anharmonicity into the bath, thus giving an effectively harmonic CL model with a frequency of (for an application to water, see e.g. Fried1992 ()). As a word of caution, we note that this model is tailored to yield linear response correctly, whereas other observables are not expected to be reproduced.
The third route is based on the mapping , which employs the non-linear projection technique. Note that although the derivation is exact, the memory kernel has been linearized to restore the desired equation structure. The resulting spectra, orange curves in Fig. 2, improve slightly with respect to that of the direct mapping (green lines) in terms of peak positions, but still have a wrong shape compared to the reference spectra (black lines). Hence, the mapping is not sufficient to completely describe the spectra for the systems studied, and the only successful approach is the effectively harmonic mapping . As a side-result we note that intramolecular forces appear to be not very important, as it follows from the comparison of panels a) and b) in Fig. 2.
To summarize, we have considered three ways to construct a mapping and have shown that only , based on an effectively harmonic potential, gives correct spectra. However, such an effective harmonic system potential would not be suitable for describing nonlinear infrared spectroscopy since nonlinear response functions will vanish in this case. In principle, there could be a yet unknown mapping yielding anharmonic system potentials. To this end we would like to consider the problem from a more general perspective without being bound to a particular way of constructing the mapping.
We note that so far our assessment was based on vibrational spectra, that is on . In fact, performing the dynamics according to Eq. (2) yields two presumably independent correlation functions, and thus defines a mapping . In order to have a self-consistent formalism the pair of correlation functions ; should coincide with the pair ; obtained from the explicit MD simulations that was used to establish mapping . Mathematically, we hence require to be the inverse mapping of . Unfortunately, mapping is not invertible, since it constitutes a mapping of two functions onto one, unless the two functions are mutually dependent. Thus, cannot be the inverse mapping of in the general case, and the self-consistency requirement of the formalism is generally corrupted. This problem is termed the “invertibility problem”.
To illustrate this point on a particular example, the procedure described above led to two strikingly different pairs of for mapping , as is shown in Fig. 3 for the case of HOD in water. Importantly, in the case of the mapping this problem does not surface, since for a harmonic force and are mutually dependent as . Thus, the mapping can be formulated with alone and the invertibility problem does not occur. This suggests that it is the anharmonicity of the system potential, , in the GLE that is at the root of all deficiencies of the applicability of the CL model to real anharmonic molecular systems. Since this consideration is independent on the particular form of the mapping , the invertibility problem would surface in general, that is, for an arbitrary anharmonic system.
Conclusions and Outlook
In this Letter the applicability of the CL model to the dynamics of realistic anharmonic systems has been investigated. Three possible mappings between a real system and the model, i.e. the direct mapping and the mappings based on (non-)linear projection techniques, have been considered. Investigating three representative systems, we obtained numerical evidence that only the mapping based on the linear projection technique () reproduces the linear absorption spectrum. Importantly, this comes at the price of projecting the system’s anharmonicity into the bath, resulting in an effectively harmonic potential which is not suitable for describing nonlinear spectroscopy. The general argument against a mapping that leads to an anharmonic system potential and gives a reasonable description have been provided and illustrated on the particular example of the direct mapping. This argument, termed invertibility problem, is not bound to the particular way of constructing the mapping and thus can be considered as a strong indication that such a mapping does not exist within the (bi-linear) CL model.
One might ask whether extending the model by considering a non-linear coupling on the system side would yield a solution. Unfortunately, the structure of the GLE derived from such an extended model would differ from that obtained from the rigorous projection techniques, see Supplement for a particular example of square-linear model. This makes establishing the mapping for an extended model difficult if at all possible.
The authors gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (Sfb 652). We thank Dr. Sergey Bokarev for fruitful discussions.
- (1) V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, 3rd Revised and Enlarged Edition, Wiley-VCH, Weinheim, 2011.
- (2) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and M. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
- (3) A. O. Caldeira and A. J. Leggett, Phys. Rev. Lett. 49, 1545 (1982).
- (4) A. O. Caldeira and A. J. Leggett, Ann. Phys. (N. Y). 149, 374 (1983).
- (5) H. Grabert, P. Schramm, and G.-L. Ingold, Phys. Rep. 168, 115 (1988).
- (6) S. Mukamel, Principles of Nonlinear Optical Spectroscopy, Oxford University Press, 1995.
- (7) Y. Tanimura and S. Mukamel, Phys. Rev. E 47, 118 (1993).
- (8) R. Zwanzig, Nonequilibrium Statistical Mechanics, Oxford University Press, 2001.
- (9) G. W. Ford and M. Kac, J. Stat. Phys. 46, 803 (1987).
- (10) G. W. Ford, J. T. Lewis, and R. F. Oconnell, Phys. Rev. A 37, 4419 (1988).
- (11) R. Feynman and F. Vernon, Ann. Phys. (N. Y). 24, 118 (1963).
- (12) A. Caldeira and A. Leggett, Phys. A Stat. Mech. its Appl. 121, 587 (1983).
- (13) A. G. Dijkstra and Y. Tanimura, Phys. Rev. Lett. 104, 250401 (2010).
- (14) D. Suess, A. Eisfeld, and W. T. Strunz, Phys. Rev. Lett. 113, 150403 (2014).
- (15) K. Giese, M. Petković, H. Naundorf, and O. Kühn, Phys. Rep. 430, 211 (2006).
- (16) S. A. Egorov, E. Rabani, and B. J. Berne, J. Chem. Phys. 110, 5238 (1999).
- (17) M. Ceriotti, G. Bussi, and M. Parrinello, Phys. Rev. Lett. 102, 020601 (2009).
- (18) M. Ceriotti, G. Bussi, and M. Parrinello, J. Chem. Theory Comput. 6, 1170 (2010).
- (19) O. Kühn and Y. Tanimura, J. Chem. Phys. 119, 2155 (2003).
- (20) Y. Tanimura and A. Ishizaki, Acc. Chem. Res. 42, 1270 (2009).
- (21) N. Makri, J. Phys. Chem. B 103, 2823 (1999).
- (22) G. Goodyear and R. M. Stratt, J. Chem. Phys. 105, 10050 (1996).
- (23) M. Tuckerman and B. J. Berne, J. Chem. Phys. 98, 7301 (1993).
- (24) H. Mori, Prog. Theor. Phys. 33, 423 (1965).
- (25) D. Cubero and S. N. Yaliraki, J. Chem. Phys. 122, 034108 (2005).
- (26) S. Kawai and T. Komatsuzaki, J. Chem. Phys. 134, 114523 (2011).
- (27) O. F. Lange and H. Grubmüller, J. Chem. Phys. 124, 214903 (2006).
- (28) J. Stenger, D. Madsen, P. Hamm, E. T. J. Nibbering, and T. Elsaesser, Phys. Rev. Lett. 87, 027401 (2001).
- (29) C. Roth, S. Chatzipapadopoulos, D. Kerlé, F. Friedriszik, M. Lütgens, S. Lochbrunner, O. Kühn, and R. Ludwig, New J. Phys. 14, 105026 (2012).
- (30) F. Paesani and G. A. Voth, J. Chem. Phys. 132, 014105 (2010).
- (31) T. Köddermann, D. Paschek, and R. Ludwig, ChemPhysChem 8, 2464 (2007).
- (32) B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008).
- (33) S. D. Ivanov, A. Witt, and D. Marx, Phys. Chem. Chem. Phys. 15, 10270 (2013).
- (34) L. E. Fried, N. Bernstein, and S. Mukamel, Phys. Rev. Lett. 68, 1842 (1992).