Quantum effects in intermediate-temperature dipole-dipole correlation-functions in the presence of an environment F. Grossmann Citation: The Journal of Chemical Physics 141, 144305 (2014); doi: 10.1063/1.4896835 View online: http://dx.doi.org/10.1063/1.4896835 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/14?ver=pdfcov Published by the AIP Publishing Articles you may be interested in MULTIMODE quantum calculations of vibrational energies and IR spectrum of the NO+(H2O) cluster using accurate potential energy and dipole moment surfaces J. Chem. Phys. 141, 124311 (2014); 10.1063/1.4896200 Molecular near-field antenna effect in resonance hyper-Raman scattering: Intermolecular vibronic intensity borrowing of solvent from solute through dipole-dipole and dipole-quadrupole interactions J. Chem. Phys. 140, 204506 (2014); 10.1063/1.4879058 Quantum theory of atoms in molecules/charge-charge flux-dipole flux models for fundamental vibrational intensity changes on H-bond formation of water and hydrogen fluoride J. Chem. Phys. 140, 084306 (2014); 10.1063/1.4865938 Effects of permanent dipole moments in transient four-wave mixing experiments J. Chem. Phys. 127, 094107 (2007); 10.1063/1.2753472 A variational study of nuclear dynamics and structural flexibility of the CH 2 OH radical J. Chem. Phys. 119, 3098 (2003); 10.1063/1.1591730
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
THE JOURNAL OF CHEMICAL PHYSICS 141, 144305 (2014)
Quantum effects in intermediate-temperature dipole-dipole correlation-functions in the presence of an environment F. Grossmanna) Institut für Theoretische Physik, Technische Universität Dresden, 01062 Dresden, Germany
(Received 25 August 2014; accepted 19 September 2014; published online 9 October 2014) We study thermal dipole-dipole correlation functions and their corresponding IR spectra in the presence of an intermediate temperature bath of harmonic oscillators. Whereas for a Morse oscillator without coupling to a heat bath, the quantum level structure is displayed in the spectrum, classical calculations show a broad, smeared out spectrum. In the presence of purely Ohmic dissipation already for extremely weak relaxation rate, the classical case is recovered. Using the HEOM approach of Tanimura and Wolynes [Phys. Rev. A 43, 4131 (1991)], we show that to observe some remnants of the level structure in the spectrum for moderate damping strength as well as intermediate temperatures, the dynamics has to be non-Markovian (Ohmic dissipation with cutoff). © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4896835] I. INTRODUCTION
Quantum mechanical correlation functions play a key role in the theoretical understanding of the interaction between matter and radiation or particles. Such diverse phenomena as the absorption of light by molecules and the inelastic scattering of neutrons from solids are described by correlation functions and their Fourier transforms. A recent review of that field, including both linear and nonlinear response theories is given in Ref. 1. In the present work, we focus our investigation on molecular vibrational spectroscopy in the presence of a thermal heat bath. The theoretical description of vibrational dephasing has been in the focus of interest from an early stage2 and models with a stochastic modulation of the vibrational energies have been used,3, 4 that have critical drawbacks because thermal equilibrium is never reached.1 A theoretical understanding of infrared molecular spectra relies on the evaluation of dipole-dipole correlation functions.5 Approaches using imaginary time path integral techniques6 suffer from instabilities of analytic continuation,7 while molecular dynamics methods8, 9 fail to capture quantum effects. Linearized semiclassical methods10 as well as centroid molecular dynamics (CMD)11 and the ring polymer molecular dynamics (RPMD) method12 have been used recently to calculate dipole-dipole correlation functions for complex molecular systems, like liquid water. Investigations of the applicability of RPMD12 as well as both RPMD and centroid molecular dynamics for IR spectra13 have been given. A solution to the contamination problem of RPMD by internal modes of the ring polymer has been addressed recently.14 Here, the influence of a thermal bath of harmonic oscillators on the dynamics of an infrared-active system is taken into account by implementing the corresponding partial differential equation analog of the Feynman-Vernon functional integral description of the reduced density dynamics.15 For a) Electronic mail: [email protected]
0021-9606/2014/141(14)/144305/7/$30.00
the general (non-Markovian) case with Ohmic dissipation and a Lorentzian cutoff, the corresponding reduced hierarchy of equations of motion (HEOM) is based on earlier work by Tanimura and Kubo16 and, for bilinear system bath coupling, has been given in the form used herein by Tanimura and Wolynes.1, 17 In the high temperature Markovian (white noise) limit their approach simplifies to the (single) master equation of Caldeira and Leggett.18 An analogous hierarchy of coupled density matrix equations of motion has been introduced by Meier and Tannor.19 It uses a parametrization of the spectral density of bath oscillators and is not restricted to a special form of the system bath interaction. In standard applications of IR spectroscopy, the anharmonicity of the bond does not play a central role, due to the fact that at room temperature or below, mostly only the vibrational ground state of the system studied is populated and a single line will be observed in the spectrum. Peaks from the population of higher lying vibrational levels are also present but strongly suppressed.20 In recent experiments on molecules at surfaces, however, also higher vibrational states are appreciably populated and in so-called sum frequency generation spectroscopy also the red shifted peaks coming from the higher lying states are observed.21 Very recently, research has directly focused on the IR spectroscopy (of CO2 ) at high temperatures, up to 6000 K.22, 23 By using high temperature thermal density matrices as initial states, we will directly calculate the corresponding dipole-dipole correlation function in different levels of description: full quantum and linearized semiclassical without coupling to the bath and HEOM and Caldeira-Leggett master equation in the presence of the coupling. The paper is organized as follows: A brief introduction to correlation functions for IR spectroscopy is given in Sec. II. The numerical methods used to evaluate the dipole-dipole correlation function are sketched in Sec. III. In Sec. IV, we show and discuss numerical results in the cases without and with system bath coupling. We conclude the presentation by hinting at possible future developments.
141, 144305-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
144305-2
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
II. CORRELATION FUNCTIONS FOR IR SPECTRA
A molecular system, weakly probed by an external light source can be treated using Fermi’s golden rule of perturbation theory. Using the Heisenberg picture of quantum mechanics, the (IR)-absorption line shape for an isotropic absorbing system can be expressed as5, 24 ∞ 1 I (ω) = dt C(t) exp{−iωt}, (1) 2π −∞
III. NUMERICAL METHODS
with the dipole-dipole correlation function C(t) = μ · μ(t),
In this expression, a mid-point rule for the argument of the potential has been used. An end point rule or an average potential rule may lead to appreciable differences in error estimates for propagators27 but the results would be barely different for any of the numerics below. The good quality of the high temperature approximation also in the case of intermediate temperatures can be seen in Sec. IV by comparison to the stick spectra results in the case without damping.
(2)
and the total electric dipole moment operator μ of the system under study (operator hats suppressed here and in the following). The brackets indicate thermal equilibrium averaging. The real-valuedness of the spectrum is guaranteed by the fact that the real (imaginary) part of the correlation function is an even (odd) function of time
In most of the numerics, we choose to use cyclic permutation under the trace in (5), to let all the dynamical factors act onto the dipole-modified density matrix in the fashion 1 i i μμ(t) = dxx|μe− ¯ H t e−βH μe ¯ H t |x Z 1 = dxx|μρ(x, ˜ x, t)|x (8) Z
Re C(t) = Re C(−t),
(3)
with ρ˜ = ρ μ.
Im C(t) = −Im C(−t).
(4)
A. Numerics without damping
In the position representation, and in 1D for a single IR active mode, the correlation function becomes 1 i i μμ(t) = (5) dxx|μe ¯ H t μe− ¯ H t e−βH |x Z with the partition function Z and β = 1/kT being proportional to the inverse of the temperature T. For the 1D harmonic oscillator, V (x) = 12 mω02 x 2 , one finds the exact analytical result ¯ω0 β ¯ cos(ω0 t) coth + i sin(ω0 t) . μμ(t) = 2mω0 2 (6) Both the real as well as the imaginary part oscillate with the same frequency. The spectrum therefore consists of a single line at ω = ω0 . For more realistic, anharmonic potentials one generally does not have an exact analytic expression for the correlation function at hand and we will use different numerical approaches to its calculation below. From the absorption line shape, the imaginary part of the dielectric constant can be calculated by multiplication with a (known) frequency dependent factor.25 The thermal density matrix for the calculation of the equilibrium average in (5) can be gained from an imaginary time propagation. We do not want to complicate matters, however, but would like to have an analytical initial state at hand. This is possible by restriction to the case of high temperatures. Furthermore, it will turn out that only for high temperatures, interesting quantum level structure will be visible in the IR spectrum. For high temperatures, i.e., for small β, a short time approximation to the “imaginary time propagator” exp (−βH) is possible. It leads to Ref. 26 x+x 0 m 1 − m 2 (x−x0 )2 −βV 2 2πβ¯ e e . (7) ρβ (x, x0 ) ≈ 2 Z 2πβ¯
1. Full quantum case
The numerical method employed to calculate the propagated ρ˜ in the quantum case without damping is the split operator method,28 using “thermal wavepackets” at half β in the following way: First we split the density operator into two parts at twice the temperature and insert unity in terms of position eigenstates 3 times into (5). The density matrix elements β
β (x0 , x) ≡ x0 |e− 2 H |x, β
˜ β (x, x0 ) ≡ x|e− 2 H μ|x0 ,
(9) (10)
that appear in the equation that emerges can be considered as “wavefunctions” and we think of the argument x as a parameter to be integrated over. Due to the imaginary time splitting of the Boltzmann operator into two parts, the high temperature approximation becomes more accurate. We note that in the numerics with damping (see below), we use the (unsplitted) density matrix and if that would be used also for the calculation of spectra in the undamped case only minute differences of the peaks at low frequencies can be observed (not shown). The unitary time-evolution operators in (5) propagate these thermal wavefunctions forward, respectively backward in time and we thus arrive at the expression 1 ˜ β∗ (x0 , x, t)μ(x0 )β (x0 , x, t) dx dx0 μμ(t) = Z (11) for the correlation function. We neglect the factor 1/Z in the numerics to be presented below, because it just leads to the multiplication of the spectrum with an overall factor. 2. The classical case
In the purely classical case, for the undamped results, we employed the linearized semiclassical initial value
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
144305-3
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
representation.29 This method is also known as the classical Wigner method and we sketch its application to the calculation of dipole-dipole correlation functions briefly. The classical Wigner method is based on the linearized semiclassical IVR (LSC-IVR) expression 1 β dqdpAW (q, p)BW (qt , pt ) (12) C LSC−IVR (t) = 2π ¯ that can be gained from the double phase space integral that emerges after inserting two semiclassical IVR expressions for the propagators forward and backward in time in (5) and a subsequent linearization of the exponent.29 ˆ ˆ the calFor the operators Aˆ = Z −1 e−β H qˆ and Bˆ = q, culation of the Wigner transforms is trivial in the hightemperature case and leads to the expression Zcl−1 LSC−IVR C (t) = (13) dqdpe−βH (p,q) q q(t) 2π ¯ for the dipole-dipole correlation function (we neglect the (classical) partition function in the numerics). Here, the dynamics is still acting on the dipole and the equilibrium density is not propagated. In passing, we note that both forward-backward semiclassical methods, the one developed by the Makri group30 as well as the (non-generalized) one of the Miller group31, 32 in the present case of dipole-dipole correlation functions do not go beyond the LSC-IVR level of description, i.e., they are not able to account for interference effects.33 B. Numerics with damping
It is commonly believed that the presence of a heat bath enforces a transition to classicality in the dynamics of the damped system, see, e.g., Ref. 34. In the present case, this transition can be studied by coupling the IR-active mode to a Caldeira-Leggett bath of harmonic oscillators, leading to the Hamiltonian18 H = HS + HB + HSB ,
(14)
which is a sum of three terms, representing the infrared active system, a heat bath (reservoir), and the interaction of the system with the reservoir, which we assume here to be bilinear. Other cases like the square-linear one have also been studied in a similar context.35, 36 Apart from its temperature, the decisive quantity specifying the heat bath is its spectral density. Frequently, a linear spectral density with Lorentzian cutoff of the form37 J (ω) = mγ ω
2 ωD 2 ωD + ω2
(15)
is used. Here γ is the relaxation rate (or damping strength) and ωD is the so-called Drude frequency. Feynman and Vernon15 studied the (reduced) system dynamics for the Hamiltonian in (14) in the path integral language for the case of a factorized initial condition (i.e., the initial density matrix of the composite system is taken as a product of a system density matrix times a bath density matrix). They identified the effect of the environment on the re-
duced density matrix of the system in terms of a so-called influence functional. Caldeira and Leggett (CL)18 have later derived a corresponding master equation (a partial differential equation (PDE) analog of the path integral) for the reduced density matrix elements if the spectral density of the oscillators of the heat bath is Ohmic, i.e., strictly linear in frequency (ωD → ∞ in (15)) and in the high temperature limit. The dynamics of the reduced density matrix of the infrared active system (here assumed to be one-dimensional) is then Markovian (memory free) and follows the master equation (we set m and ¯ equal to unity in the following)
γ ρ(x, ˙ x , t) = −iLA − (x − x )(∂x − ∂x ) 2 γ (16) − (x − x )2 ρ(x, x , t) β with 1 1 LA = − ∂x2 + ∂x2 + V (x) − V (x ), 2 2
(17)
where the second and third term on the RHS of (16) are the dissipation and the thermal fluctuation term, respectively. In a phase space representation, this equation has later on also been used for intermediate temperatures, with surprisingly good agreement with numerically exact path integral results,38 as well as for the description of two-dimensional spectroscopy.39 We have implemented it in position space, using an implicit alternating direction scheme.40, 41 In the numerics we are using the analytical high-temperature density matrix (7) modified with the dipole operator but without the normalizing factor 1/Z as the initial state for the propagation. In addition, we have also implemented the KleinKramers42, 43 Fokker-Planck equation ρ˙W (p, q, t) = {HS , ρW }Poi + γ ∂p [pρW (p, q, t)] γ + ∂p2 ρW (p, q, t) β
(18)
in Wigner phase space, where {. . . }Poi indicates the Poisson bracket. This is the classical analog of the CL equation and for zero damping (γ = 0) gives the same result as the LSC-IVR from above. The assumption of a linearly increasing Ohmic density of bath oscillators is rather unphysical (leading to a white noise driving in a corresponding Langevin equation), however. We have therefore, in addition, implemented an approach going beyond the Markovian master equation by Caldeira and Leggett. This more general approach was devised by Tanimura and Wolynes,17 based on earlier work by Tanimura and Kubo,16 to describe in principle arbitrary temperature baths (for low temperatures, care has to be taken, see the Appendix of Ref. 17) with the spectral density given in (15). An analogous generalization was later also given by Meier and Tannor for parametrized spectral densities.19 To account for memory effects, the PDE for the reduced density matrix, now denoted by ρ 0 , has to be augmented by auxiliary density matrices, ρ n , n = 1, 2,. . . , leading to a hierarchy of coupled PDEs (reduced
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
0.6 0.4 0.2 0 -0.2 -0.4
0.03
(a) I(ω)
Re C(t)
144305-4
0
50
100 t
150
0.02 0.01 0
200
(b)
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 1. The real part of the full quantum dipole-dipole correlation (a) and IR line shape (b) of a 1D Morse oscillator (both in arbitrary units) for dimensionless temperature T = 5. Sticks in (b) follow from (23), continuous line from FT of full quantum dipole-dipole correlation function; the continuous line curve has been scaled by an overall factor such that the absolute maximum coincides with that of the stick spectrum.
hierarchy equations of motion: HEOM)
IV. NUMERICAL RESULTS FOR THE IR LINE SHAPE
ρ˙n (x, x , t) = −[iLA +nωD ]ρn (x, x , t)−i(x −x )ρn+1 (x, x , t) βωD γ ωD
(∂x −∂x ) + ωD cot (x −x ) −ni 2 2
A. The one-dimensional IR-active model system
× ρn−1 (x, x , t).
(19)
This system of equations is truncated in a numerical application by a so-called terminator equation βωD γ
ρN+1 (x, x , t) ≈ −i (∂x −∂x ) + ωD cot (x −x ) 2 2 × ρN (x, x , t),
(20)
closing the hierarchy of N + 1 recurrence differential equations. For large ωD , one may set N = 0 in (19) and for T → ∞, it has been shown that the approach is equivalent to the one by Caldeira and Leggett.17 Following Tanimura and Wolynes, we have implemented the HEOM using an explicit fourth order Runge Kutta (RK) method for the time step and finite difference expressions for the spatial derivatives.
0.6 0.4 0.2 0 -0.2 -0.4
V (x) = D[1 − exp(−αx)]2 .
0.03
(a)
0
20
40
60 t
80
100
(21)
The dimensionless potential parameters that we use are √ D = 100 and α = 0.2 2, leading to ω0 = 4 for the frequency of (harmonic) oscillations around the minimum and a total of 50 bound states, typical for diatomic molecules.45 In contrast to the harmonic oscillator, the ladder of eigenenergies is not equidistantly spaced but the spacings decrease with increasing energy of the Morse oscillator. The effect this has for the IR-spectrum is seen in Fig. 1. There we compare two spectra at an intermediate temperature of T = 5, corresponding to around 3000 K in a typical molecular system (CO-stretch). We used this intermediate temperature, because it allows for an appreciable population of excited vibrational levels. Even higher temperatures would lead to classical results already for very low damping strengths, however.
I(ω)
C(t)
The model system we study is a 1D Morse oscillator44 with unit mass and the potential
(b)
0.02 0.01 0
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 2. The classical dipole correlation function (a) and the IR line shape (b) of a 1D Morse oscillator (both in arbitrary units) for dimensionless temperature T = 5. Sticks in (b) follow from (23), continuous line from FT of LSC-IVR dipole-dipole correlation function; the continuous curve has been scaled with an overall factor such that it touches the largest stick’s maximum.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
144305-5
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
I(ω)
0.16 0.12 0.08 0.04 0
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 3. The IR line-shape of a 1D Morse oscillator (in arbitrary units) for dimensionless temperature T = 2. Sticks follow from (23), continuous line from FT of LSC-IVR dipole-dipole correlation function and multiplied with an overall factor to match the largest stick.
The sticks in Fig. 1 are calculated by using energy eigenstates for the evaluation of the correlation function, 1 i|μ|j j |μ|iei(Ej −Ei )t e−βEi , (22) C(t) = Z i,j which upon Fourier transformation gives 1 I (ω) = |i|μ|j |2 δ[ω − (Ej − Ei )]e−βEi Z i,j
(23)
0.6 0.4 0.2 0 -0.2 -0.4
B. Morse oscillator coupled to a bath of harmonic oscillators
In the dissipative case, we solved the Tanimura-Wolynes HEOM for nonzero as well as for zero N, corresponding to the Caldeira-Leggett case. The number of grid points we used are 512 × 512 (although almost identical results can be gained by using a 400 × 400 grid) for a grid interval of x ∈ [−5, 15] and the time step in the RK scheme was chosen as 0.001. Implementing the CL case with an implicit scheme, this time step could be chosen much larger (typically 0.005). We start with the time series for the dipole-dipole correlation in the CL case and again for an intermediate value of T = 5 for two different damping strengths, displayed in panel (a) of Fig. 4. We do not display the corresponding classical results of the Klein-Kramers Fokker Planck equation, because they are almost identical to the CL ones for the present parameters (a similar observation has been made in Ref. 41). In contrast to panel (a) of Fig. 1 but very similarly to panel (a) of Fig. 2, both results quickly decay to zero after t = 50, there is no recurrence of the oscillation amplitude. The corresponding IR line-shapes are shown in panel (b) of Figure 4. No trace of the quantum stick spectrum is visible
0.02
(a) I(ω)
Re C(t)
for the line shape. We then restrict the calculation of the double sum to terms with j = i + 1 and use the Morse oscillator energies44 (from which we can also calculate the partition function) and the dipole matrix elements given by Gallas.46 The finite width spectra can be derived from the Fourier-transform (FT) of the correlation function, calculated numerically with split operator FFT. The finite width of those Fourier transform results (performed with zero padding and windowing47 ) has no physical reason in the present case but is due to the finite length of the corresponding time signal. The very good agreement of the stick spectrum calculated without any approximation and the spectrum following from the thermal wavefunctions that uses a high temperature approximation to the initial state underpins the high quality of
this approximation for the present parameters (see also the remark above about the imaginary time splitting of the initial state). A similar comparison with the stick spectrum, but now with the spectrum gained by FT of the classical LSC-IVR signal, shown in panel (a) of Fig. 2 is given in panel (b) of Fig. 2. As to be expected, the classical dynamics does not “know” about the quantum level structure of the Morse oscillator but just shows the decrease of the oscillation frequency with increasing energy in a smooth fashion. The envelope over the stick spectrum is reproduced by the classical result. Furthermore, we checked this fact also for lower temperatures and found that the lack of zero-point energy in the classical result is reflected in the fact that high frequency contributions (coming from low energy oscillations around the minimum of the oscillator) are overestimated severely at low T as can be seen in Fig. 3.
0
20
40
60 t
80
100
(b)
0.01
0
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 4. (a) Real part of the dipole-dipole correlation (in arbitrary units) of a 1D Morse oscillator for dimensionless temperature T = 5 coupled to an Ohmic bath with γ = 0.05 (solid black line), γ = 0.1 (dashed green line). (b) Corresponding line shapes (in arbitrary units).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
0.6 0.4 0.2 0 -0.2 -0.4
0.04
(a)
(b)
0.03 I(ω)
Re C(t)
144305-6
0.02 0.01
0
20
40
60
80
100
t
0
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 5. (a) Real part of the dipole-dipole correlation (in arbitrary units) of a 1D Morse oscillator for dimensionless temperature T = 5 coupled to an Ohmic bath with Lorentzian cutoff (ωD = 1) with γ = 0.05 (solid black line), γ = 0.1 (dashed green line). (b) Corresponding line shapes (arbitrary units).
any more, also the continuous curve does not act as the envelope of the (undamped) stick spectrum any more. Increasing the dissipation strength does, however, red-shift the maximum of the line shape. This effect, however, is not a quantum effect but is also present in the data derived from the Klein-Kramers equation (not shown). The complete washing out of quantum effects in the CL results (no sub-peaks in the line-shape any more) sets in already at very small damping strengths, around γ = 0.01 (not shown). In order to still see some quantum effects in the presence of intermediate coupling strength, the bath has to be more realistic, i.e., the high frequency components of the spectral density have to be suppressed. Using the full HEOM, we therefore now study the case of a Drude-Lorentz cutoff of the spectral density and choose the Drude frequency to be below the system frequency by setting ωD = 1, making resonant energy transfer between system and bath unlikely by effectively suppressing the spectral density at resonance by an order of magnitude. This choice of parameter is typical for systems of chemical interest20 and has also been used in a semiclassical study of interference quenching, using Ohmic damping with an exponential cutoff.48 Furthermore, we have used N = 5 for the termination of the hierarchy and have checked that higher N does not change the results. In Fig. 5, the results for two different damping strengths are shown. In contrast to the undamped classical as well as the CL case, the dynamics in panel (a) shows a recurrence of the amplitude of the correlation function. For the spectra in panel (b) of Fig. 5, in strong contrast to the pure Ohmic dissipation case, now, even for intermediate γ = 0.1, a remnant to the stick spectrum in Fig. 1 is still visible. V. CONCLUSIONS AND OUTLOOK
Quantum effects are prone to vanish as soon as the (quantum) system of interest is coupled to an environment. We have explicitly studied the case of the line shape of infrared absorption in the presence of a bath of intermediate temperature. We stress that the CL master equation, strictly, is not applicable to this case, due to the fact that it was derived for high temperatures (T → ∞) where quantum effects do go away already for extremely small coupling strength (due to the decoher-
ing influence of the fluctuation term proportional to γ /β). It can, however, be applied for finite temperatures as well. Even in this case, already for intermediate damping strength, due to the purely Ohmic spectral density, quantum effects do get washed out. The work of Tanimura and co-workers has opened a road to consistently treat also intermediate temperatures, as well as more realistic spectral densities of the bath oscillators. Even for the case of low temperatures, Sakurai and Tanimura have observed a small 1-2 peak in the IR spectrum of a typical vibrational mode of a peptide.20 Here, we have shown that for intermediate temperature, as well as intermediate damping strengths, to observe quantum effects in the line shape (remnants of stick spectra), the effective damping strength at resonance has to be reduced (Drude-Lorentz cutoff) and the full-fledged HEOM has to be considered, allowing for the non-Markovian description of the dynamics of the reduced density matrix. The necessity of a non-Markovian approach has also been stressed for the correct treatment of nonlinear vibrational line shapes by Yang et al.49 Recently, work with HEOM has focused on the identification of quantum effects in two dimensional spectroscopy1, 20, 50–52 as well as in (photosynthetic) electronic energy transfer.53, 54 With the IR line-shape calculations shown here, it becomes clear that for intermediate temperatures and anharmonic potentials even in this much simpler type of experiment some quantum effects may still be seen if the damping strength can be kept low. In future work, it would be highly desirable to understand the effect that initial correlation between the system and the bath would have on the spectra. A first account of the treatment of initial correlations in the HEOM context has been given recently.55 Furthermore, also a generalization to more general coupling schemes, e.g., square-linear,1, 20 instead of bilinear system bath coupling or more general spectral densities of the bath oscillators are worthwhile topics for future investigations. ACKNOWLEDGMENTS
F.G. gratefully acknowledges valuable discussions with Professor Y. Tanimura, Professor A. Ishizaki, and Dr. N.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
144305-7
F. Grossmann
Dattani. Furthermore, financial support from the Japan Society for the Promotion of Science in the form of a shortterm fellowship to visit Kyoto University is also gratefully acknowledged. 1 Y.
Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006). G. Rothschild, J. Chem. Phys. 65, 455 (1976). 3 D. Oxtoby, Adv. Chem. Phys. 40, 1 (1979). 4 D. Oxtoby, Adv. Chem. Phys. 47, 487 (1981). 5 R. G. Gordon, Adv. Magn. Res. 3, 1 (1968). 6 D. Thirumalai and B. J. Berne, J. Chem. Phys. 79, 5029 (1983). 7 E. Gallicchio, S. A. Egorov, and B. J. Berne, J. Chem. Phys. 109, 7745 (1998). 8 P. H. Behrens and K. R. Wilson, J. Chem. Phys. 74, 4872 (1981). 9 B. M. Auer and J. L. Skinner, J. Chem. Phys. 128, 224511 (2008). 10 J. Liu, W. H. Miller, F. Paesani, W. Zhang, and D. A. Case, J. Chem. Phys. 131, 164509 (2009); J. Liu, W. H. Miller, G. S. Fanourgakis, S. S. Xantheas, S. Imoto, and S. Saito, J. Chem. Phys. 135, 244503 (2011). 11 F. Paesani, S. S. Xantheas, and G. A. Voth, J. Phys. Chem. B 113, 13118 (2009). 12 S. Habershon and D. E. Manolopoulos, J. Chem. Phys. 131, 244518 (2009). 13 A. Witt, S. D. Ivanov, M. Shiga, H. Forbert, and D. Marx, J. Chem. Phys. 130, 194510 (2009). 14 M. Rossi, M. Ceriotti, and D. E. Manolopouolos, J. Chem. Phys. 140, 234116 (2014). 15 R. P. Feynman and F. L. Vernon, Ann. Phys. (N.Y.) 24, 118 (1963). 16 Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn. 58, 101 (1989). 17 Y. Tanimura and P. G. Wolynes, Phys. Rev. A 43, 4131 (1991). 18 A. O. Caldeira and A. J. Leggett, Physica A 121, 587 (1983). 19 C. Meier and D. Tannor, J. Chem. Phys. 111, 3365 (1999). 20 A. Sakurai and Y. Tanimura, J. Phys. Chem. A 115, 4009 (2011). 21 M. Bonn, C. Hess, and M. Wolf, J. Chem. Phys. 115, 7725 (2001). 22 S. Depraz, M. Y. Perrin, and A. Soufiani, J. Quant. Spectrosc. Radiat. Transfer 113, 1 (2012). 23 S. Depraz, M. Y. Perrin, Ph. Rivière, and A. Soufiani, J. Quant. Spectrosc. Radiat. Transfer 113, 14 (2012). 24 E. J. Heller, Acc. Chem. Res. 14, 368 (1981). 25 R. Ramrez, T. López-Ciudad, P. Kumar P, and D. Marx, J. Chem. Phys. 121, 3973 (2004). 2 W.
J. Chem. Phys. 141, 144305 (2014) 26 R.
P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals (Mc Graw-Hill, New York, 1965). 27 N. Makri and W. H. Miller, Chem. Phys. Lett. 151, 1 (1988). 28 M. D. Feit, J. A. Fleck Jr., and A. Steiger, J. Comput. Phys. 47, 412 (1982). 29 H. Wang, X. Sun, and W. H. Miller, J. Chem. Phys. 108, 9726 (1998). 30 K. Thomson and N. Makri, J. Chem. Phys. 110, 1343 (1999). 31 W. H. Miller, Faraday Discuss. Chem. Soc. 110, 1 (1998). 32 X. Sun and W. H. Miller, J. Chem. Phys. 110, 6635 (1999). 33 M. Thoss, H. Wang, and W. H. Miller, J. Chem. Phys. 114, 9220 (2001). 34 M. Kryvohuz and J. Cao, J. Chem. Phys. 130, 234107 (2009), and references therein. 35 T. Kato and Y. Tanimura, J. Chem. Phys. 117, 6221 (2002). 36 A. Ishizaki and Y. Tanimura, J. Chem. Phys. 125, 084501 (2006). 37 U. Weiss, Quantum Dissipative Systems, 2nd ed. (World Scientific, 1999). 38 S. Zhang and E. Pollak, J. Chem. Phys. 118, 4357 (2003). 39 Y. Tanimura, Chem. Phys. 233, 217 (1998). 40 F. Grossmann, J. Chem. Phys. 103, 3696 (1995). 41 F. Grossmann and W. Koch, J. Chem. Phys. 130, 034105 (2009). 42 O. Klein, Ark. Mat. Astron. Fys. 16, 1 (1921). 43 H. A. Kramers, Physica 7, 284 (1940). 44 P. M. Morse, Phys. Rev. 34, 57 (1929). 45 G. Karl, Phys. Rev. A 49, 799 (1994), note that in order to be infrared active, diatomic molecules must be unsymmetric (e.g., heteronuclear). 46 J. A. C. Gallas, Phys. Rev. A 21, 1829 (1980). 47 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran, 2nd ed. (Cambridge University Press, Cambridge, 1992). 48 H. Wang, M. Thoss, K. L. Sorge, R. Gelabert, X. Giménez, and W. H. Miller, J. Chem. Phys. 114, 2562 (2001). 49 S. Yang, J. Shao, and J. Cao, J. Chem. Phys. 121, 11250 (2004). 50 S. Mukamel, Y. Tanimura, and P. Hamm, Acc. Chem. Res. 42, 1207 (2009). 51 L. Chen, R. Zheng, Q. Shi, and Y. Yan, J. Chem. Phys. 132, 024505 (2010). 52 B. Hein, C. Kreisbeck, T. Kramer, and M. Rodríguez, New J. Phys. 14, 023018 (2012). 53 A. Ishizaki and G. R. Fleming, Proc. Natl. Acad. Sci. U.S.A. 106, 17255 (2009). 54 Y. Tanimura, J. Chem. Phys. 137, 22A550 (2012). 55 Y. Tanimura, J. Chem. Phys. 141, 044114 (2014).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
THE JOURNAL OF CHEMICAL PHYSICS 141, 144305 (2014)
Quantum effects in intermediate-temperature dipole-dipole correlation-functions in the presence of an environment F. Grossmanna) Institut für Theoretische Physik, Technische Universität Dresden, 01062 Dresden, Germany
(Received 25 August 2014; accepted 19 September 2014; published online 9 October 2014) We study thermal dipole-dipole correlation functions and their corresponding IR spectra in the presence of an intermediate temperature bath of harmonic oscillators. Whereas for a Morse oscillator without coupling to a heat bath, the quantum level structure is displayed in the spectrum, classical calculations show a broad, smeared out spectrum. In the presence of purely Ohmic dissipation already for extremely weak relaxation rate, the classical case is recovered. Using the HEOM approach of Tanimura and Wolynes [Phys. Rev. A 43, 4131 (1991)], we show that to observe some remnants of the level structure in the spectrum for moderate damping strength as well as intermediate temperatures, the dynamics has to be non-Markovian (Ohmic dissipation with cutoff). © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4896835] I. INTRODUCTION
Quantum mechanical correlation functions play a key role in the theoretical understanding of the interaction between matter and radiation or particles. Such diverse phenomena as the absorption of light by molecules and the inelastic scattering of neutrons from solids are described by correlation functions and their Fourier transforms. A recent review of that field, including both linear and nonlinear response theories is given in Ref. 1. In the present work, we focus our investigation on molecular vibrational spectroscopy in the presence of a thermal heat bath. The theoretical description of vibrational dephasing has been in the focus of interest from an early stage2 and models with a stochastic modulation of the vibrational energies have been used,3, 4 that have critical drawbacks because thermal equilibrium is never reached.1 A theoretical understanding of infrared molecular spectra relies on the evaluation of dipole-dipole correlation functions.5 Approaches using imaginary time path integral techniques6 suffer from instabilities of analytic continuation,7 while molecular dynamics methods8, 9 fail to capture quantum effects. Linearized semiclassical methods10 as well as centroid molecular dynamics (CMD)11 and the ring polymer molecular dynamics (RPMD) method12 have been used recently to calculate dipole-dipole correlation functions for complex molecular systems, like liquid water. Investigations of the applicability of RPMD12 as well as both RPMD and centroid molecular dynamics for IR spectra13 have been given. A solution to the contamination problem of RPMD by internal modes of the ring polymer has been addressed recently.14 Here, the influence of a thermal bath of harmonic oscillators on the dynamics of an infrared-active system is taken into account by implementing the corresponding partial differential equation analog of the Feynman-Vernon functional integral description of the reduced density dynamics.15 For a) Electronic mail: [email protected]
0021-9606/2014/141(14)/144305/7/$30.00
the general (non-Markovian) case with Ohmic dissipation and a Lorentzian cutoff, the corresponding reduced hierarchy of equations of motion (HEOM) is based on earlier work by Tanimura and Kubo16 and, for bilinear system bath coupling, has been given in the form used herein by Tanimura and Wolynes.1, 17 In the high temperature Markovian (white noise) limit their approach simplifies to the (single) master equation of Caldeira and Leggett.18 An analogous hierarchy of coupled density matrix equations of motion has been introduced by Meier and Tannor.19 It uses a parametrization of the spectral density of bath oscillators and is not restricted to a special form of the system bath interaction. In standard applications of IR spectroscopy, the anharmonicity of the bond does not play a central role, due to the fact that at room temperature or below, mostly only the vibrational ground state of the system studied is populated and a single line will be observed in the spectrum. Peaks from the population of higher lying vibrational levels are also present but strongly suppressed.20 In recent experiments on molecules at surfaces, however, also higher vibrational states are appreciably populated and in so-called sum frequency generation spectroscopy also the red shifted peaks coming from the higher lying states are observed.21 Very recently, research has directly focused on the IR spectroscopy (of CO2 ) at high temperatures, up to 6000 K.22, 23 By using high temperature thermal density matrices as initial states, we will directly calculate the corresponding dipole-dipole correlation function in different levels of description: full quantum and linearized semiclassical without coupling to the bath and HEOM and Caldeira-Leggett master equation in the presence of the coupling. The paper is organized as follows: A brief introduction to correlation functions for IR spectroscopy is given in Sec. II. The numerical methods used to evaluate the dipole-dipole correlation function are sketched in Sec. III. In Sec. IV, we show and discuss numerical results in the cases without and with system bath coupling. We conclude the presentation by hinting at possible future developments.
141, 144305-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
144305-2
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
II. CORRELATION FUNCTIONS FOR IR SPECTRA
A molecular system, weakly probed by an external light source can be treated using Fermi’s golden rule of perturbation theory. Using the Heisenberg picture of quantum mechanics, the (IR)-absorption line shape for an isotropic absorbing system can be expressed as5, 24 ∞ 1 I (ω) = dt C(t) exp{−iωt}, (1) 2π −∞
III. NUMERICAL METHODS
with the dipole-dipole correlation function C(t) = μ · μ(t),
In this expression, a mid-point rule for the argument of the potential has been used. An end point rule or an average potential rule may lead to appreciable differences in error estimates for propagators27 but the results would be barely different for any of the numerics below. The good quality of the high temperature approximation also in the case of intermediate temperatures can be seen in Sec. IV by comparison to the stick spectra results in the case without damping.
(2)
and the total electric dipole moment operator μ of the system under study (operator hats suppressed here and in the following). The brackets indicate thermal equilibrium averaging. The real-valuedness of the spectrum is guaranteed by the fact that the real (imaginary) part of the correlation function is an even (odd) function of time
In most of the numerics, we choose to use cyclic permutation under the trace in (5), to let all the dynamical factors act onto the dipole-modified density matrix in the fashion 1 i i μμ(t) = dxx|μe− ¯ H t e−βH μe ¯ H t |x Z 1 = dxx|μρ(x, ˜ x, t)|x (8) Z
Re C(t) = Re C(−t),
(3)
with ρ˜ = ρ μ.
Im C(t) = −Im C(−t).
(4)
A. Numerics without damping
In the position representation, and in 1D for a single IR active mode, the correlation function becomes 1 i i μμ(t) = (5) dxx|μe ¯ H t μe− ¯ H t e−βH |x Z with the partition function Z and β = 1/kT being proportional to the inverse of the temperature T. For the 1D harmonic oscillator, V (x) = 12 mω02 x 2 , one finds the exact analytical result ¯ω0 β ¯ cos(ω0 t) coth + i sin(ω0 t) . μμ(t) = 2mω0 2 (6) Both the real as well as the imaginary part oscillate with the same frequency. The spectrum therefore consists of a single line at ω = ω0 . For more realistic, anharmonic potentials one generally does not have an exact analytic expression for the correlation function at hand and we will use different numerical approaches to its calculation below. From the absorption line shape, the imaginary part of the dielectric constant can be calculated by multiplication with a (known) frequency dependent factor.25 The thermal density matrix for the calculation of the equilibrium average in (5) can be gained from an imaginary time propagation. We do not want to complicate matters, however, but would like to have an analytical initial state at hand. This is possible by restriction to the case of high temperatures. Furthermore, it will turn out that only for high temperatures, interesting quantum level structure will be visible in the IR spectrum. For high temperatures, i.e., for small β, a short time approximation to the “imaginary time propagator” exp (−βH) is possible. It leads to Ref. 26 x+x 0 m 1 − m 2 (x−x0 )2 −βV 2 2πβ¯ e e . (7) ρβ (x, x0 ) ≈ 2 Z 2πβ¯
1. Full quantum case
The numerical method employed to calculate the propagated ρ˜ in the quantum case without damping is the split operator method,28 using “thermal wavepackets” at half β in the following way: First we split the density operator into two parts at twice the temperature and insert unity in terms of position eigenstates 3 times into (5). The density matrix elements β
β (x0 , x) ≡ x0 |e− 2 H |x, β
˜ β (x, x0 ) ≡ x|e− 2 H μ|x0 ,
(9) (10)
that appear in the equation that emerges can be considered as “wavefunctions” and we think of the argument x as a parameter to be integrated over. Due to the imaginary time splitting of the Boltzmann operator into two parts, the high temperature approximation becomes more accurate. We note that in the numerics with damping (see below), we use the (unsplitted) density matrix and if that would be used also for the calculation of spectra in the undamped case only minute differences of the peaks at low frequencies can be observed (not shown). The unitary time-evolution operators in (5) propagate these thermal wavefunctions forward, respectively backward in time and we thus arrive at the expression 1 ˜ β∗ (x0 , x, t)μ(x0 )β (x0 , x, t) dx dx0 μμ(t) = Z (11) for the correlation function. We neglect the factor 1/Z in the numerics to be presented below, because it just leads to the multiplication of the spectrum with an overall factor. 2. The classical case
In the purely classical case, for the undamped results, we employed the linearized semiclassical initial value
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
144305-3
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
representation.29 This method is also known as the classical Wigner method and we sketch its application to the calculation of dipole-dipole correlation functions briefly. The classical Wigner method is based on the linearized semiclassical IVR (LSC-IVR) expression 1 β dqdpAW (q, p)BW (qt , pt ) (12) C LSC−IVR (t) = 2π ¯ that can be gained from the double phase space integral that emerges after inserting two semiclassical IVR expressions for the propagators forward and backward in time in (5) and a subsequent linearization of the exponent.29 ˆ ˆ the calFor the operators Aˆ = Z −1 e−β H qˆ and Bˆ = q, culation of the Wigner transforms is trivial in the hightemperature case and leads to the expression Zcl−1 LSC−IVR C (t) = (13) dqdpe−βH (p,q) q q(t) 2π ¯ for the dipole-dipole correlation function (we neglect the (classical) partition function in the numerics). Here, the dynamics is still acting on the dipole and the equilibrium density is not propagated. In passing, we note that both forward-backward semiclassical methods, the one developed by the Makri group30 as well as the (non-generalized) one of the Miller group31, 32 in the present case of dipole-dipole correlation functions do not go beyond the LSC-IVR level of description, i.e., they are not able to account for interference effects.33 B. Numerics with damping
It is commonly believed that the presence of a heat bath enforces a transition to classicality in the dynamics of the damped system, see, e.g., Ref. 34. In the present case, this transition can be studied by coupling the IR-active mode to a Caldeira-Leggett bath of harmonic oscillators, leading to the Hamiltonian18 H = HS + HB + HSB ,
(14)
which is a sum of three terms, representing the infrared active system, a heat bath (reservoir), and the interaction of the system with the reservoir, which we assume here to be bilinear. Other cases like the square-linear one have also been studied in a similar context.35, 36 Apart from its temperature, the decisive quantity specifying the heat bath is its spectral density. Frequently, a linear spectral density with Lorentzian cutoff of the form37 J (ω) = mγ ω
2 ωD 2 ωD + ω2
(15)
is used. Here γ is the relaxation rate (or damping strength) and ωD is the so-called Drude frequency. Feynman and Vernon15 studied the (reduced) system dynamics for the Hamiltonian in (14) in the path integral language for the case of a factorized initial condition (i.e., the initial density matrix of the composite system is taken as a product of a system density matrix times a bath density matrix). They identified the effect of the environment on the re-
duced density matrix of the system in terms of a so-called influence functional. Caldeira and Leggett (CL)18 have later derived a corresponding master equation (a partial differential equation (PDE) analog of the path integral) for the reduced density matrix elements if the spectral density of the oscillators of the heat bath is Ohmic, i.e., strictly linear in frequency (ωD → ∞ in (15)) and in the high temperature limit. The dynamics of the reduced density matrix of the infrared active system (here assumed to be one-dimensional) is then Markovian (memory free) and follows the master equation (we set m and ¯ equal to unity in the following)
γ ρ(x, ˙ x , t) = −iLA − (x − x )(∂x − ∂x ) 2 γ (16) − (x − x )2 ρ(x, x , t) β with 1 1 LA = − ∂x2 + ∂x2 + V (x) − V (x ), 2 2
(17)
where the second and third term on the RHS of (16) are the dissipation and the thermal fluctuation term, respectively. In a phase space representation, this equation has later on also been used for intermediate temperatures, with surprisingly good agreement with numerically exact path integral results,38 as well as for the description of two-dimensional spectroscopy.39 We have implemented it in position space, using an implicit alternating direction scheme.40, 41 In the numerics we are using the analytical high-temperature density matrix (7) modified with the dipole operator but without the normalizing factor 1/Z as the initial state for the propagation. In addition, we have also implemented the KleinKramers42, 43 Fokker-Planck equation ρ˙W (p, q, t) = {HS , ρW }Poi + γ ∂p [pρW (p, q, t)] γ + ∂p2 ρW (p, q, t) β
(18)
in Wigner phase space, where {. . . }Poi indicates the Poisson bracket. This is the classical analog of the CL equation and for zero damping (γ = 0) gives the same result as the LSC-IVR from above. The assumption of a linearly increasing Ohmic density of bath oscillators is rather unphysical (leading to a white noise driving in a corresponding Langevin equation), however. We have therefore, in addition, implemented an approach going beyond the Markovian master equation by Caldeira and Leggett. This more general approach was devised by Tanimura and Wolynes,17 based on earlier work by Tanimura and Kubo,16 to describe in principle arbitrary temperature baths (for low temperatures, care has to be taken, see the Appendix of Ref. 17) with the spectral density given in (15). An analogous generalization was later also given by Meier and Tannor for parametrized spectral densities.19 To account for memory effects, the PDE for the reduced density matrix, now denoted by ρ 0 , has to be augmented by auxiliary density matrices, ρ n , n = 1, 2,. . . , leading to a hierarchy of coupled PDEs (reduced
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
0.6 0.4 0.2 0 -0.2 -0.4
0.03
(a) I(ω)
Re C(t)
144305-4
0
50
100 t
150
0.02 0.01 0
200
(b)
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 1. The real part of the full quantum dipole-dipole correlation (a) and IR line shape (b) of a 1D Morse oscillator (both in arbitrary units) for dimensionless temperature T = 5. Sticks in (b) follow from (23), continuous line from FT of full quantum dipole-dipole correlation function; the continuous line curve has been scaled by an overall factor such that the absolute maximum coincides with that of the stick spectrum.
hierarchy equations of motion: HEOM)
IV. NUMERICAL RESULTS FOR THE IR LINE SHAPE
ρ˙n (x, x , t) = −[iLA +nωD ]ρn (x, x , t)−i(x −x )ρn+1 (x, x , t) βωD γ ωD
(∂x −∂x ) + ωD cot (x −x ) −ni 2 2
A. The one-dimensional IR-active model system
× ρn−1 (x, x , t).
(19)
This system of equations is truncated in a numerical application by a so-called terminator equation βωD γ
ρN+1 (x, x , t) ≈ −i (∂x −∂x ) + ωD cot (x −x ) 2 2 × ρN (x, x , t),
(20)
closing the hierarchy of N + 1 recurrence differential equations. For large ωD , one may set N = 0 in (19) and for T → ∞, it has been shown that the approach is equivalent to the one by Caldeira and Leggett.17 Following Tanimura and Wolynes, we have implemented the HEOM using an explicit fourth order Runge Kutta (RK) method for the time step and finite difference expressions for the spatial derivatives.
0.6 0.4 0.2 0 -0.2 -0.4
V (x) = D[1 − exp(−αx)]2 .
0.03
(a)
0
20
40
60 t
80
100
(21)
The dimensionless potential parameters that we use are √ D = 100 and α = 0.2 2, leading to ω0 = 4 for the frequency of (harmonic) oscillations around the minimum and a total of 50 bound states, typical for diatomic molecules.45 In contrast to the harmonic oscillator, the ladder of eigenenergies is not equidistantly spaced but the spacings decrease with increasing energy of the Morse oscillator. The effect this has for the IR-spectrum is seen in Fig. 1. There we compare two spectra at an intermediate temperature of T = 5, corresponding to around 3000 K in a typical molecular system (CO-stretch). We used this intermediate temperature, because it allows for an appreciable population of excited vibrational levels. Even higher temperatures would lead to classical results already for very low damping strengths, however.
I(ω)
C(t)
The model system we study is a 1D Morse oscillator44 with unit mass and the potential
(b)
0.02 0.01 0
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 2. The classical dipole correlation function (a) and the IR line shape (b) of a 1D Morse oscillator (both in arbitrary units) for dimensionless temperature T = 5. Sticks in (b) follow from (23), continuous line from FT of LSC-IVR dipole-dipole correlation function; the continuous curve has been scaled with an overall factor such that it touches the largest stick’s maximum.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
144305-5
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
I(ω)
0.16 0.12 0.08 0.04 0
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 3. The IR line-shape of a 1D Morse oscillator (in arbitrary units) for dimensionless temperature T = 2. Sticks follow from (23), continuous line from FT of LSC-IVR dipole-dipole correlation function and multiplied with an overall factor to match the largest stick.
The sticks in Fig. 1 are calculated by using energy eigenstates for the evaluation of the correlation function, 1 i|μ|j j |μ|iei(Ej −Ei )t e−βEi , (22) C(t) = Z i,j which upon Fourier transformation gives 1 I (ω) = |i|μ|j |2 δ[ω − (Ej − Ei )]e−βEi Z i,j
(23)
0.6 0.4 0.2 0 -0.2 -0.4
B. Morse oscillator coupled to a bath of harmonic oscillators
In the dissipative case, we solved the Tanimura-Wolynes HEOM for nonzero as well as for zero N, corresponding to the Caldeira-Leggett case. The number of grid points we used are 512 × 512 (although almost identical results can be gained by using a 400 × 400 grid) for a grid interval of x ∈ [−5, 15] and the time step in the RK scheme was chosen as 0.001. Implementing the CL case with an implicit scheme, this time step could be chosen much larger (typically 0.005). We start with the time series for the dipole-dipole correlation in the CL case and again for an intermediate value of T = 5 for two different damping strengths, displayed in panel (a) of Fig. 4. We do not display the corresponding classical results of the Klein-Kramers Fokker Planck equation, because they are almost identical to the CL ones for the present parameters (a similar observation has been made in Ref. 41). In contrast to panel (a) of Fig. 1 but very similarly to panel (a) of Fig. 2, both results quickly decay to zero after t = 50, there is no recurrence of the oscillation amplitude. The corresponding IR line-shapes are shown in panel (b) of Figure 4. No trace of the quantum stick spectrum is visible
0.02
(a) I(ω)
Re C(t)
for the line shape. We then restrict the calculation of the double sum to terms with j = i + 1 and use the Morse oscillator energies44 (from which we can also calculate the partition function) and the dipole matrix elements given by Gallas.46 The finite width spectra can be derived from the Fourier-transform (FT) of the correlation function, calculated numerically with split operator FFT. The finite width of those Fourier transform results (performed with zero padding and windowing47 ) has no physical reason in the present case but is due to the finite length of the corresponding time signal. The very good agreement of the stick spectrum calculated without any approximation and the spectrum following from the thermal wavefunctions that uses a high temperature approximation to the initial state underpins the high quality of
this approximation for the present parameters (see also the remark above about the imaginary time splitting of the initial state). A similar comparison with the stick spectrum, but now with the spectrum gained by FT of the classical LSC-IVR signal, shown in panel (a) of Fig. 2 is given in panel (b) of Fig. 2. As to be expected, the classical dynamics does not “know” about the quantum level structure of the Morse oscillator but just shows the decrease of the oscillation frequency with increasing energy in a smooth fashion. The envelope over the stick spectrum is reproduced by the classical result. Furthermore, we checked this fact also for lower temperatures and found that the lack of zero-point energy in the classical result is reflected in the fact that high frequency contributions (coming from low energy oscillations around the minimum of the oscillator) are overestimated severely at low T as can be seen in Fig. 3.
0
20
40
60 t
80
100
(b)
0.01
0
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 4. (a) Real part of the dipole-dipole correlation (in arbitrary units) of a 1D Morse oscillator for dimensionless temperature T = 5 coupled to an Ohmic bath with γ = 0.05 (solid black line), γ = 0.1 (dashed green line). (b) Corresponding line shapes (in arbitrary units).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
F. Grossmann
J. Chem. Phys. 141, 144305 (2014)
0.6 0.4 0.2 0 -0.2 -0.4
0.04
(a)
(b)
0.03 I(ω)
Re C(t)
144305-6
0.02 0.01
0
20
40
60
80
100
t
0
3
3.2 3.4 3.6 3.8 ω
4
4.2
FIG. 5. (a) Real part of the dipole-dipole correlation (in arbitrary units) of a 1D Morse oscillator for dimensionless temperature T = 5 coupled to an Ohmic bath with Lorentzian cutoff (ωD = 1) with γ = 0.05 (solid black line), γ = 0.1 (dashed green line). (b) Corresponding line shapes (arbitrary units).
any more, also the continuous curve does not act as the envelope of the (undamped) stick spectrum any more. Increasing the dissipation strength does, however, red-shift the maximum of the line shape. This effect, however, is not a quantum effect but is also present in the data derived from the Klein-Kramers equation (not shown). The complete washing out of quantum effects in the CL results (no sub-peaks in the line-shape any more) sets in already at very small damping strengths, around γ = 0.01 (not shown). In order to still see some quantum effects in the presence of intermediate coupling strength, the bath has to be more realistic, i.e., the high frequency components of the spectral density have to be suppressed. Using the full HEOM, we therefore now study the case of a Drude-Lorentz cutoff of the spectral density and choose the Drude frequency to be below the system frequency by setting ωD = 1, making resonant energy transfer between system and bath unlikely by effectively suppressing the spectral density at resonance by an order of magnitude. This choice of parameter is typical for systems of chemical interest20 and has also been used in a semiclassical study of interference quenching, using Ohmic damping with an exponential cutoff.48 Furthermore, we have used N = 5 for the termination of the hierarchy and have checked that higher N does not change the results. In Fig. 5, the results for two different damping strengths are shown. In contrast to the undamped classical as well as the CL case, the dynamics in panel (a) shows a recurrence of the amplitude of the correlation function. For the spectra in panel (b) of Fig. 5, in strong contrast to the pure Ohmic dissipation case, now, even for intermediate γ = 0.1, a remnant to the stick spectrum in Fig. 1 is still visible. V. CONCLUSIONS AND OUTLOOK
Quantum effects are prone to vanish as soon as the (quantum) system of interest is coupled to an environment. We have explicitly studied the case of the line shape of infrared absorption in the presence of a bath of intermediate temperature. We stress that the CL master equation, strictly, is not applicable to this case, due to the fact that it was derived for high temperatures (T → ∞) where quantum effects do go away already for extremely small coupling strength (due to the decoher-
ing influence of the fluctuation term proportional to γ /β). It can, however, be applied for finite temperatures as well. Even in this case, already for intermediate damping strength, due to the purely Ohmic spectral density, quantum effects do get washed out. The work of Tanimura and co-workers has opened a road to consistently treat also intermediate temperatures, as well as more realistic spectral densities of the bath oscillators. Even for the case of low temperatures, Sakurai and Tanimura have observed a small 1-2 peak in the IR spectrum of a typical vibrational mode of a peptide.20 Here, we have shown that for intermediate temperature, as well as intermediate damping strengths, to observe quantum effects in the line shape (remnants of stick spectra), the effective damping strength at resonance has to be reduced (Drude-Lorentz cutoff) and the full-fledged HEOM has to be considered, allowing for the non-Markovian description of the dynamics of the reduced density matrix. The necessity of a non-Markovian approach has also been stressed for the correct treatment of nonlinear vibrational line shapes by Yang et al.49 Recently, work with HEOM has focused on the identification of quantum effects in two dimensional spectroscopy1, 20, 50–52 as well as in (photosynthetic) electronic energy transfer.53, 54 With the IR line-shape calculations shown here, it becomes clear that for intermediate temperatures and anharmonic potentials even in this much simpler type of experiment some quantum effects may still be seen if the damping strength can be kept low. In future work, it would be highly desirable to understand the effect that initial correlation between the system and the bath would have on the spectra. A first account of the treatment of initial correlations in the HEOM context has been given recently.55 Furthermore, also a generalization to more general coupling schemes, e.g., square-linear,1, 20 instead of bilinear system bath coupling or more general spectral densities of the bath oscillators are worthwhile topics for future investigations. ACKNOWLEDGMENTS
F.G. gratefully acknowledges valuable discussions with Professor Y. Tanimura, Professor A. Ishizaki, and Dr. N.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
144305-7
F. Grossmann
Dattani. Furthermore, financial support from the Japan Society for the Promotion of Science in the form of a shortterm fellowship to visit Kyoto University is also gratefully acknowledged. 1 Y.
Tanimura, J. Phys. Soc. Jpn. 75, 082001 (2006). G. Rothschild, J. Chem. Phys. 65, 455 (1976). 3 D. Oxtoby, Adv. Chem. Phys. 40, 1 (1979). 4 D. Oxtoby, Adv. Chem. Phys. 47, 487 (1981). 5 R. G. Gordon, Adv. Magn. Res. 3, 1 (1968). 6 D. Thirumalai and B. J. Berne, J. Chem. Phys. 79, 5029 (1983). 7 E. Gallicchio, S. A. Egorov, and B. J. Berne, J. Chem. Phys. 109, 7745 (1998). 8 P. H. Behrens and K. R. Wilson, J. Chem. Phys. 74, 4872 (1981). 9 B. M. Auer and J. L. Skinner, J. Chem. Phys. 128, 224511 (2008). 10 J. Liu, W. H. Miller, F. Paesani, W. Zhang, and D. A. Case, J. Chem. Phys. 131, 164509 (2009); J. Liu, W. H. Miller, G. S. Fanourgakis, S. S. Xantheas, S. Imoto, and S. Saito, J. Chem. Phys. 135, 244503 (2011). 11 F. Paesani, S. S. Xantheas, and G. A. Voth, J. Phys. Chem. B 113, 13118 (2009). 12 S. Habershon and D. E. Manolopoulos, J. Chem. Phys. 131, 244518 (2009). 13 A. Witt, S. D. Ivanov, M. Shiga, H. Forbert, and D. Marx, J. Chem. Phys. 130, 194510 (2009). 14 M. Rossi, M. Ceriotti, and D. E. Manolopouolos, J. Chem. Phys. 140, 234116 (2014). 15 R. P. Feynman and F. L. Vernon, Ann. Phys. (N.Y.) 24, 118 (1963). 16 Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn. 58, 101 (1989). 17 Y. Tanimura and P. G. Wolynes, Phys. Rev. A 43, 4131 (1991). 18 A. O. Caldeira and A. J. Leggett, Physica A 121, 587 (1983). 19 C. Meier and D. Tannor, J. Chem. Phys. 111, 3365 (1999). 20 A. Sakurai and Y. Tanimura, J. Phys. Chem. A 115, 4009 (2011). 21 M. Bonn, C. Hess, and M. Wolf, J. Chem. Phys. 115, 7725 (2001). 22 S. Depraz, M. Y. Perrin, and A. Soufiani, J. Quant. Spectrosc. Radiat. Transfer 113, 1 (2012). 23 S. Depraz, M. Y. Perrin, Ph. Rivière, and A. Soufiani, J. Quant. Spectrosc. Radiat. Transfer 113, 14 (2012). 24 E. J. Heller, Acc. Chem. Res. 14, 368 (1981). 25 R. Ramrez, T. López-Ciudad, P. Kumar P, and D. Marx, J. Chem. Phys. 121, 3973 (2004). 2 W.
J. Chem. Phys. 141, 144305 (2014) 26 R.
P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals (Mc Graw-Hill, New York, 1965). 27 N. Makri and W. H. Miller, Chem. Phys. Lett. 151, 1 (1988). 28 M. D. Feit, J. A. Fleck Jr., and A. Steiger, J. Comput. Phys. 47, 412 (1982). 29 H. Wang, X. Sun, and W. H. Miller, J. Chem. Phys. 108, 9726 (1998). 30 K. Thomson and N. Makri, J. Chem. Phys. 110, 1343 (1999). 31 W. H. Miller, Faraday Discuss. Chem. Soc. 110, 1 (1998). 32 X. Sun and W. H. Miller, J. Chem. Phys. 110, 6635 (1999). 33 M. Thoss, H. Wang, and W. H. Miller, J. Chem. Phys. 114, 9220 (2001). 34 M. Kryvohuz and J. Cao, J. Chem. Phys. 130, 234107 (2009), and references therein. 35 T. Kato and Y. Tanimura, J. Chem. Phys. 117, 6221 (2002). 36 A. Ishizaki and Y. Tanimura, J. Chem. Phys. 125, 084501 (2006). 37 U. Weiss, Quantum Dissipative Systems, 2nd ed. (World Scientific, 1999). 38 S. Zhang and E. Pollak, J. Chem. Phys. 118, 4357 (2003). 39 Y. Tanimura, Chem. Phys. 233, 217 (1998). 40 F. Grossmann, J. Chem. Phys. 103, 3696 (1995). 41 F. Grossmann and W. Koch, J. Chem. Phys. 130, 034105 (2009). 42 O. Klein, Ark. Mat. Astron. Fys. 16, 1 (1921). 43 H. A. Kramers, Physica 7, 284 (1940). 44 P. M. Morse, Phys. Rev. 34, 57 (1929). 45 G. Karl, Phys. Rev. A 49, 799 (1994), note that in order to be infrared active, diatomic molecules must be unsymmetric (e.g., heteronuclear). 46 J. A. C. Gallas, Phys. Rev. A 21, 1829 (1980). 47 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran, 2nd ed. (Cambridge University Press, Cambridge, 1992). 48 H. Wang, M. Thoss, K. L. Sorge, R. Gelabert, X. Giménez, and W. H. Miller, J. Chem. Phys. 114, 2562 (2001). 49 S. Yang, J. Shao, and J. Cao, J. Chem. Phys. 121, 11250 (2004). 50 S. Mukamel, Y. Tanimura, and P. Hamm, Acc. Chem. Res. 42, 1207 (2009). 51 L. Chen, R. Zheng, Q. Shi, and Y. Yan, J. Chem. Phys. 132, 024505 (2010). 52 B. Hein, C. Kreisbeck, T. Kramer, and M. Rodríguez, New J. Phys. 14, 023018 (2012). 53 A. Ishizaki and G. R. Fleming, Proc. Natl. Acad. Sci. U.S.A. 106, 17255 (2009). 54 Y. Tanimura, J. Chem. Phys. 137, 22A550 (2012). 55 Y. Tanimura, J. Chem. Phys. 141, 044114 (2014).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.174.255.116 On: Mon, 01 Dec 2014 16:53:25
