The modular decomposition of the path integral is a linear-scaling, numerically exact algorithm for calculating dynamical properties of extended systems composed of multilevel units with local couplings. In a recent article, we generalized the method to wavefunction propagation in aggregates characterized by non-diagonal couplings between adjacent units. Here, we extend the method to the calculation of reduced density matrices in aggregates where each unit includes an arbitrary number of coupled harmonic bath modes, which may describe intramolecular normal mode vibrations, at finite temperature. The effects of harmonic modes are included through influence functional factors, which involve analytical expressions that we derive. Representative applications to spin arrays described by the Heisenberg Hamiltonian with dissipative interactions and to J-aggregates of perylene bisimide, where all coupled normal modes are treated explicitly, are presented.

For several decades, Feynman’s path integral approach to time-dependent quantum mechanics1,2 has found wide application in all areas of theoretical physics. In combination with Monte Carlo sampling methods,3 its imaginary formulation4 has been used as a numerically exact and efficient simulation tool for computing equilibrium properties of quantum many-body systems obeying Boltzmann or Bose statistics.5,6 Such simulations have led to a wealth of information regarding the structure of low-temperature fluids, superfluids, and the Bose–Einstein condensate.7 Furthermore, the quantum–classical isomorphism8 motivated the introduction of molecular dynamics techniques to equilibrium path integral calculations,9–12 which are more efficient for some types of interactions.

The main obstacle in applying Feynman’s formulation to real-time dynamics is the presence of a highly oscillatory phase, which arises from the complex exponential of the action integral along each quantum path.13 The paradigm of quantum dissipation is based on the interaction of a system with a harmonic bath,14 where a large (often infinite) number of harmonic degrees of freedom can mimic many effects of generic condensed phase environments. For example, the effects of dissipation on tunneling dynamics15 and reaction rate constants16,17 have been thoroughly investigated this way. The effective harmonic bath picture has also been shown to arise as the dominant contribution from a condensed phase environment in the so-called Gaussian (or linear) response limit18 and justifies the parabolic free energy surface model employed in electron transfer theory.19 Recently, rigorous quantum–classical path integral20 (QCPI) simulations of the ferrocene–ferrocenium electron transfer in a solvent with strongly nonlinear intermolecular forces21 established the quantitative accuracy of the effective harmonic bath picture for some processes involving charge transfer. Furthermore, harmonic degrees of freedom arise naturally in the description of intramolecular normal modes, and the reaction path Hamiltonian22 employs exactly this type of system–bath interaction.

A distinct advantage of the path integral approach is that harmonic bath degrees of freedom may be integrated out analytically within this formulation, giving rise to an influence functional23 that captures exactly the effects of the harmonic bath on the dynamics of the quantum system of interest. However, the influence functional is nonlocal in time, and thus, the path integral cannot be evaluated by Markovian matrix-vector procedures. The iterative quasi-adiabatic propagator path integral (i-QuAPI) is based on tensor decompositions of the reduced density matrix,24,25 two-time,26,27 or multi-time correlation functions,28 which capture bath-induced memory effects and allow for propagation to long times. The full i-QuAPI algorithm involves the storage of tensors whose size grows exponentially with the memory length, but path filtering procedures,29,30 coarse graining,31 blip,32,33 and singular value decomposition (SVD)-based matrix product operator (MPO) decompositions34 can dramatically reduce the cost. A recent small matrix disentanglement of the path integral variables35,36 (SMatPI) reduces the required storage to that of ordinary matrices of size equal to that of the bare system, extending the feasibility of numerically exact simulation to multi-state systems.

Unlike path integral-based methods, approaches based on the Schrödinger equation have to treat the bath degrees of freedom explicitly. Among wavefunction-based schemes, the multiconfiguration time-dependent Hartree (MCTDH) method37–39 along with its multilayer40 and effective mode41 extensions has shown impressive success on select system–bath models.42 The main advantage of the MCTDH approach is its ability to also treat anharmonic bath degrees of freedom. However, wavefunction-based methods are naturally suited to zero-temperature properties, and their extension to finite temperature by summing many microcanonical calculations is computationally prohibitive when the environment includes a large number of low-frequency modes. In the special case of a harmonic bath described by a spectral density of the Drude form, the hierarchical equations of motion (HEOM) method43 allows efficient simulation of system–bath dynamics. However, its use in more general situations requires a decomposition of the spectral density into Drude-type components,44 and since the method scales exponentially with respect to the number of such terms, this approach becomes impractical for simulating processes in real, structured environments. Methods based on the density matrix renormalization group45,46 (DMRG) and matrix product states47,48 can be highly efficient for extended systems, but their cost increases steeply with the number of bath degrees of freedom, and thus, these methods are restricted to a small number of vibrational modes. A recent development combining matrix product techniques with a pseudomode representation of the bath49 appears to overcome some of these limitations in cases where the bath spectral density can be efficiently decomposed, but its ability to tackle molecular environments with structured vibrational spectra may be limited.

Recent developments have extended the applicability of real-time path integral methods to more complex situations. In particular, it has been shown that extended systems described in terms of discrete Hamiltonians with sufficiently local interactions can be treated exactly via a modular path integral (MPI) approach,50,51 where the quantum paths of a unit are linked to those of the adjacent unit and summed without encountering a sign problem. The MPI methodology is applicable to systems with a one-dimensional topology, such as chains of spins or molecular aggregates, but can also account for branching or ring structure, and leads to linear scaling with system size. While the structure of the MPI algorithm is simplest in Ising-type models, where the coupling is local in the single-unit basis, we recently extended the methodology to the general case where the coupling is not diagonalizable in any single-particle basis.52 We also introduced an efficient factorization of the MPI linking algorithm that leads to exponential savings.53 

In this paper, we develop an MPI methodology suitable for propagation of density matrices, where each of the units includes coupling to a finite-temperature harmonic bath, which might represent intramolecular vibrations that couple to the electronic states of the unit or a generic dissipative environment. We show that the coupled harmonic degrees of freedom can be integrated out analytically, giving rise to influence functional factors that apply separately to each unit, and derive expressions for the coefficients in the discretized influence functional expression.

We emphasize that the development described in this paper is not a method for treating system–bath dynamics. Even in the absence of bath degrees of freedom, calculations on extended quantum systems, such as the quantum Ising and Heisenberg spin models or molecular aggregates of two or more electronic states, are limited by the usual exponential scaling of quantum mechanics. For example, a basis set treatment of a chain containing r spin-12 units requires storage of a wavefunction with 2r elements, and thus, calculations quickly become impractical when r exceeds 10. On the other hand, the MPI decomposition is characterized by linear scaling with r. Moreover, the inclusion of harmonic bath degrees of freedom is possible with very little additional cost and is just as simple for discrete vibrational modes as it is for a model spectral density. These attributes suggest that the MPI methodology is able to address important problems not amenable to other approaches.

In Sec. II, we describe the aggregate Hamiltonian, where the discrete states of each unit are coupled to harmonic bath modes. For concreteness, we refer to the discrete states of each unit as “electronic states” and to the attached bath degrees of freedom as “vibrational” modes, with the understanding that a unit might be defined in terms of spin states, and the bath degrees of freedom might describe a generic dissipative environment. In Sec. III, we develop the path integral expression and the MPI algorithm for the density matrix and the expectation value of a target operator. In Sec. IV, we describe the procedure for integrating out the harmonic degrees of freedom and derive expressions for the influence functional coefficients. In Sec. V, we apply the methodology to examine the effects of a dissipative bath on the dynamics of Heisenberg spin chains and to simulate excitation energy transfer in J-aggregates of perylene bisimide, where all coupled normal modes of each monomer are treated explicitly. Finally, in Sec. VI, we give some concluding remarks.

Consider an array of r units for which the Hamiltonian has the form

(2.1)

where the subscript α = A, B, …, Y, Z labels the monomers. Each unit is described by the Hamiltonian

(2.2)

where

(2.3)

describes the Mα electronic states |nα⟩. Each of these states is coupled to vibrational degrees of freedom, which are described by the Hamiltonian

(2.4)

where

(2.5)

is the operator that specifies the coupling to the electronic system. Equation (2.4) is diagonal in the chosen electronic basis. Its diagonal element

(2.6)

describes the vibrational Hamiltonian associated with the electronic state |nα⟩. The states |nα⟩ may be eigenstates of Hαe, in which case Eq. (2.4) is the normal mode expansion of the Born–Oppenheimer potential energy surfaces, or may describe diabatic states coupled through nondiagonal terms in Eq. (2.3). The basis states might also describe spin degrees of freedom in magnetic materials or in model Hamiltonians, such as the transverse field Ising chain54–56 or the Heisenberg model,57 which interact with a local phonon or generic dissipative environment.

The molecular units are coupled to their nearest neighbors through the electronic coupling terms

(2.7)

This type of coupling is particularly relevant for modeling the electron-vibrational dynamics in molecular aggregates coupled through electrostatic or other non-covalent interactions. For example, the Frenkel exciton58 coupling form is often used to describe excitation energy transfer in molecular aggregates.

In the general case where the electronic coupling is not diagonalizable in a single-unit basis, we construct pairs of adjacent units by partitioning the Hamiltonian as follows:52 

(3.1)

where

(3.2)

etc. We use the letter Z to label the last unit of the chain and Y to label the next-to-last unit, regardless of the number of units. The expectation value of an operator O can be obtained from the propagated density matrix,

(3.3)

where ρ(0) is the initial density operator and the trace is with respect to all (electronic and vibrational) degrees of freedom of all units. The variables labeled nA,f+,,nZ,f+ and nA,f,,nZ,f indicate the endpoints of the forward and backward propagators.

In a recent paper,52 we described the MPI representation of the electronic propagator, which yields the time evolution of the electronic wavefunction. While the pair structure of the short-time propagators appears to double the number of path integral variables, we showed that the additional auxiliary variable can be eliminated by implementing a swapped-pair propagator sequence. Here, we adopt this structure and extend the MPI algorithm to the density matrix, where the units include vibrational degrees of freedom. We retain the path representation in terms of electronic states only, while the vibrational degrees of freedom appear in operator form and are immediately traced over.

A diagrammatic representation of the MPI linking expressions for the density matrix is shown in Fig. 1. Note that only even-numbered path integral variables remain in the representation of the first unit A, while odd-numbered variables survive in the last unit. To simplify the notation, we have relabeled the variables nA,f+,,nZ,f+ and nA,f,,nZ,f using a numbering scheme compatible with the path integral beads of each unit. For the first and last units, we use nA,f±nA,N± and nZ,f±nZ,N±. For units B, C, …, Y, the last path integral variable is nα,N±, so we use the notation nα,f±nα,N+1±.

FIG. 1.

Schematic illustration showing the linked variables for the first three pairs of the MPI expression for the reduced density matrix with N = 4. Each colored shape corresponds to a two-monomer short-time propagator. Colored circles represent the path integral variables. Black circles indicate the endpoints, which are connected by the matrix element of the probed operator. Ellipses represent the initial density operator. For clarity, only some of the variables of the backward paths are labeled.

FIG. 1.

Schematic illustration showing the linked variables for the first three pairs of the MPI expression for the reduced density matrix with N = 4. Each colored shape corresponds to a two-monomer short-time propagator. Colored circles represent the path integral variables. Black circles indicate the endpoints, which are connected by the matrix element of the probed operator. Ellipses represent the initial density operator. For clarity, only some of the variables of the backward paths are labeled.

Close modal

If ρ(0) is a product of single-unit operators (e.g., all spins of an Ising chain are initially in the “up” state, or a single molecule is electronically excited at t = 0 and the vibrational modes are equilibrated with respect to that state) and if the electronic properties of a particular molecular unit α1 are being probed, the expectation value of the observable operator O^=α=1nO^α, where O^α=Îα (the identity operator for monomer α) for all αα1 is given by the following sequence of MPI steps:

(3.4)

where TrA,v denotes the trace with respect to the vibrational degrees of freedom of unit A,

(3.5)
(3.6)

etc. Each of these steps links the paths in the electronic space of a monomer to those of the next unit. The B-to-C and C-to-D steps are performed sequentially until the next-to-last unit is reached. If the chain has an odd number of units, the last step is given by

(3.7)

and the observable is obtained through the sum

(3.8)

The propagators that appear in these expressions are illustrated in Fig. 1. We note that this diagram corresponds to one of the two mirror-image diagrams in the swapped pair MPI expression, which arise from reversing the order of the Trotter splitting. Averaging the results with respect to these two symmetric forms reduces the Trotter error, allowing a larger path integral time step. A factorization of the path linking operations53 leads to a very efficient algorithm that scales as LαlogMαLα, where Lα is the number of paths of unit α.

This modular, one-monomer-at-a-time decomposition of the path integral is key to ensuring a linear scaling of the computational cost with the number of monomers. Each element of the arrays Rα is a specific sequence of the path variables of monomer α and hence a particular pair of forward and backward paths for this unit. These paths only interact with those of monomer α + 1 with amplitudes that depend on the specifics of the operator Ĥα,α+1 and are linked to those through the expressions.

Consider the Hamiltonian for a unit of two adjacent monomers. The Hamiltonian for the pair α, α + 1 becomes

(4.1)

We partition this Hamiltonian into pure electronic and electron-vibration components. Noting that Hαev and Hα+1ev commute, the short-time propagator is given by the following product:

(4.2)

The path integral variables of unit A are separated by time steps 2Δt and are coupled only to those of unit B. Evaluating Eq. (4.2) in the electronic basis according to the swapped-pair representation of the path integral,52 the first two forward propagator elements for the AB pair are

(4.3)

Recall again that ĤAv(sA,nα) are operators in the space of the vibrational degrees of freedom of unit A, with the potential minimum at the coordinate values qα,i=cα,isα,nα/mωα,i2. Thus, the exponential operators in Eq. (4.3) describe evolution of a multidimensional forced harmonic oscillator, where the time dependence of the force is specified by the particular sequence of electronic states and the corresponding time lengths. Collecting all the short time propagators that correspond to the forward and backward evolution operators and including the vibrational component of the initial density operator ρAv(0), we see that the vibrational modes of unit A augment the electronic amplitude in the path integral by the factor

(4.4)

where Trv denotes the trace with the respect to the vibrational degrees of freedom and ZvA is the vibrational partition function for unit A.

The product of operators in Eq. (4.4) describes the time evolution of the vibrational degrees of freedom subject to forces determined by the value of the electron-vibration coupling index sA, which changes at multiples of 2Δt. This familiar form is the time-discretized Feynman–Vernon influence functional,23 augmented by a phase59 that arises from the “counterterms”14,15 in Eq. (2.6). If the bath is initially in thermal equilibrium at the temperature T = 1/kBβ, the influence functional for unit α has the form

(4.5)

where Cα(t′) is the force autocorrelation function of unit α and sα±(t) is the electronic index that specifies the vibrational state along the given path. The correlation function is given by the expression

(4.6)

It is often expressed in terms of the spectral density function,14 

(4.7)

To evaluate the influence functional arising from the vibrational modes of a particular unit, we discretize the double time integral, taking advantage of the piecewise constant form of the force. This brings the influence functional to the form of a discrete sum,

(4.8)

and

(4.9)

where the coefficients ηα,kk are to be determined.

In the case of unit A, the force changes at the end of each double time step, and thus, the discretization of the double integral in Eq. (4.8) is straightforward. For example, it is easy to see that

(4.10)

etc. Performing the time integrals using Eq. (4.6) leads to expressions for the discretized influence functional coefficients.

We now proceed to examine the vibrational influence functional for unit B. Equation (4.3) shows the sequence of forces that arises from the interaction with unit A during the first few time steps. However, because of the splitting of Eq. (3.2), half of the vibrational Hamiltonian of unit B is contained in the propagators for the pair BC. The first two of the forward time propagators are given by the expressions

(4.11)

Combining the vibrational operators that contain B variables in Eqs. (4.3) and (4.11), one sees that the time dependence of the force for unit B is given by the sequence specified in Fig. 2. Using this force sequence to discretize Eq. (4.5), we find

(4.12)

etc. After performing the time integrals, we obtain the discretized influence functional coefficients for the vibrational modes of unit B.

FIG. 2.

Force discretization for the vibrational influence functionals of various units in a chain. Top panel: force sequence for terminal unit Z. Middle panel: force sequence for non-terminal units. Bottom panel: force sequence for terminal unit A.

FIG. 2.

Force discretization for the vibrational influence functionals of various units in a chain. Top panel: force sequence for terminal unit Z. Middle panel: force sequence for non-terminal units. Bottom panel: force sequence for terminal unit A.

Close modal

Next, the variables of unit C enter also through the CD propagators. The first of these is split into electronic and vibrational components as follows:

(4.13)

Combining Eqs. (4.11) and (4.13), we observe that the discretization of the force for monomer C is identical to that of monomer B. In fact, it is not hard to see that every non-terminal monomer has the same discretization of the force. Therefore, all non-terminal monomers B, C, D, …, Y have an identical set of discretized influence functional coefficients, which thus need to be evaluated only once.

Finally, we obtain the third and final unique discretization of forces for unit Z. Because unit Z is only coupled to unit Y, the sequence of forces can be completely obtained from the YZ product of propagators. The first two such propagators are

(4.14)

The sequence of forces for unit Z is shown in Fig. 2. Using this force sequence, we obtain the influence functional coefficients for the last unit of the chain. The discretized vibrational influence functional coefficients are given in Tables I–III for all units.

TABLE I.

Influence functional coefficients ηα,kk,  k,k=2,4,,N2 and kk for unit A.

ηA,00 = ηA,NN 0dωJω2ω2coth12ωβ  1cosωΔt+isinωΔt 
ηA,N0 0dω2Jωω2sin212ωΔtcoth12ωβ  cosN1ωΔt+isinN1ωΔt 
ηA,kk 0dωJω2ω2coth12ωβ  1cos2ωΔt+isin2ωΔt 
ηA,k0 0dωJω2ω2sinωΔtsin12ωΔtcoth12ωβ  cosk12ωΔt+isink12ωΔt 
ηA,Nk 0dω2Jωω2sinωΔtsin12ωΔtcoth12ωβ  cosNk12ωΔt+isinNk12ωΔt 
ηA,kk 0dω2Jωω2sin2ωΔtcoth12ωβ  coskkωΔt+isinkkωΔt 
ηA,00 = ηA,NN 0dωJω2ω2coth12ωβ  1cosωΔt+isinωΔt 
ηA,N0 0dω2Jωω2sin212ωΔtcoth12ωβ  cosN1ωΔt+isinN1ωΔt 
ηA,kk 0dωJω2ω2coth12ωβ  1cos2ωΔt+isin2ωΔt 
ηA,k0 0dωJω2ω2sinωΔtsin12ωΔtcoth12ωβ  cosk12ωΔt+isink12ωΔt 
ηA,Nk 0dω2Jωω2sinωΔtsin12ωΔtcoth12ωβ  cosNk12ωΔt+isinNk12ωΔt 
ηA,kk 0dω2Jωω2sin2ωΔtcoth12ωβ  coskkωΔt+isinkkωΔt 
TABLE II.

Influence functional coefficients ηα,kk,  k,k=2,3,4,,N andkk for units α = B, C, …, Y.

ηα,00 = ηα,N+1,N+1 0dωJω2ω2coth12ωβ  1cos14ωΔt+isin1cos14ωΔt 
ηα,11 = ηα,NN 0dωJω2ω2coth12ωβ  1cos34ωΔt+isin34ωΔt 
ηα,10 = ηα,N+1,N 0dω2Jωω2sin18ωΔt  sin38ωΔtcoth12ωβ  cos12ωΔt+isin12ωΔt 
ηα,N+1,0 0dω2Jωω2sin218ωΔt   coth12ωβ  cosN14ωΔt+isinN14ωΔt 
ηα,N0 = ηα,N+1,1 0dω2Jωω2sin18ωΔt  sin38ωΔt   coth12ωβ  cosN34ωΔt+isinN34ωΔt 
ηα,N1 0dω2Jωω2sin238ωΔt   coth12ωβ  cosN54ωΔt+isinN54ωΔt 
ηα,kk 0dωJω2ω2  coth12ωβ  1cosωΔt+isinωΔt 
ηα,k0 0dω2Jωω2sin12ωΔt  sin18ωΔt   coth12ωβ  cosk58ωΔt+isink58ωΔt 
ηα,N+1,k 0dω2Jωω2sin12ωΔt  sin18ωΔt   coth12ωβ  cosNk+38ωΔt+isinNk+38ωΔt 
ηα,k1 0dω2Jωω2sin12ωΔt  sin18ωΔt   coth12ωβ  cosk98ωΔt+isink98ωΔt 
ηα,Nk 0dω2Jωω2sin12ωΔt  sin18ωΔt   coth12ωβ  cosNk18ωΔt+isinNk18ωΔt 
ηα,kk 0dω2Jωω2sin212ωΔt   coth12ωβ  coskkωΔt+isinkkωΔt 
ηα,00 = ηα,N+1,N+1 0dωJω2ω2coth12ωβ  1cos14ωΔt+isin1cos14ωΔt 
ηα,11 = ηα,NN 0dωJω2ω2coth12ωβ  1cos34ωΔt+isin34ωΔt 
ηα,10 = ηα,N+1,N 0dω2Jωω2sin18ωΔt  sin38ωΔtcoth12ωβ  cos12ωΔt+isin12ωΔt 
ηα,N+1,0 0dω2Jωω2sin218ωΔt   coth12ωβ  cosN14ωΔt+isinN14ωΔt 
ηα,N0 = ηα,N+1,1 0dω2Jωω2sin18ωΔt  sin38ωΔt   coth12ωβ  cosN34ωΔt+isinN34ωΔt 
ηα,N1 0dω2Jωω2sin238ωΔt   coth12ωβ  cosN54ωΔt+isinN54ωΔt 
ηα,kk 0dωJω2ω2  coth12ωβ  1cosωΔt+isinωΔt 
ηα,k0 0dω2Jωω2sin12ωΔt  sin18ωΔt   coth12ωβ  cosk58ωΔt+isink58ωΔt 
ηα,N+1,k 0dω2Jωω2sin12ωΔt  sin18ωΔt   coth12ωβ  cosNk+38ωΔt+isinNk+38ωΔt 
ηα,k1 0dω2Jωω2sin12ωΔt  sin18ωΔt   coth12ωβ  cosk98ωΔt+isink98ωΔt 
ηα,Nk 0dω2Jωω2sin12ωΔt  sin18ωΔt   coth12ωβ  cosNk18ωΔt+isinNk18ωΔt 
ηα,kk 0dω2Jωω2sin212ωΔt   coth12ωβ  coskkωΔt+isinkkωΔt 
TABLE III.

Influence functional coefficients ηα,kk,  k,k=2,3,4,,N andkk for unit Z.

ηZ,00 = ηZ,NN 0dωJω2ω2  coth12ωβ  1cos12ωΔt+isin12ωΔt 
ηZ,11 = ηZ,N−1,N−1 0dωJω2ω2  coth12ωβ  1cos32ωΔt+isin32ωΔt 
ηZ,10 = ηZ,N,N−1 0dωJω2ω2sin14ωΔt  sin34ωΔt  coth12ωβ  cosωΔt+isinωΔt 
ηZ,N,0 0dω2Jωω2sin214ωΔt   coth12ωβ  cosN12ωΔt+isinN12ωΔt 
ηZ,N,1 = ηZ,N−1,0 0dω2Jωω2sin14ωΔt  sin34ωΔt   coth12ωβ  cosN32ωΔt+isinN32ωΔt 
ηZ,N−1,1 0dω2Jωω2sin234ωΔt   coth12ωβ  cosN52ωΔt+isinN52ωΔt 
ηZ,kk 0dωJω2ω2  coth12ωβ  1cos2ωΔt+isin2ωΔt 
ηZ,k0 0dω2Jωω2sinωΔt  sin14ωΔt   coth12ωβ  cosk14ωΔt+isink14ωΔt 
ηZ,N,k 0dω2Jωω2sinωΔt  sin14ωΔt   coth12ωβ  cosNk14ωΔt+isinNk14ωΔt 
ηZ,k1 0dω2Jωω2sinωΔt  sin14ωΔt   coth12ωβ  cosk54ωΔt+isink54ωΔt 
ηZ,N−1,k 0dω2Jωω2sinωΔt  sin14ωΔt   coth12ωβ  cosNk54ωΔt+isinNk54ωΔt 
ηZ,kk 0dω2Jωω2sin2ωΔt     coth12ωβ  cos2kkωΔt+isin2kkωΔt 
ηZ,00 = ηZ,NN 0dωJω2ω2  coth12ωβ  1cos12ωΔt+isin12ωΔt 
ηZ,11 = ηZ,N−1,N−1 0dωJω2ω2  coth12ωβ  1cos32ωΔt+isin32ωΔt 
ηZ,10 = ηZ,N,N−1 0dωJω2ω2sin14ωΔt  sin34ωΔt  coth12ωβ  cosωΔt+isinωΔt 
ηZ,N,0 0dω2Jωω2sin214ωΔt   coth12ωβ  cosN12ωΔt+isinN12ωΔt 
ηZ,N,1 = ηZ,N−1,0 0dω2Jωω2sin14ωΔt  sin34ωΔt   coth12ωβ  cosN32ωΔt+isinN32ωΔt 
ηZ,N−1,1 0dω2Jωω2sin234ωΔt   coth12ωβ  cosN52ωΔt+isinN52ωΔt 
ηZ,kk 0dωJω2ω2  coth12ωβ  1cos2ωΔt+isin2ωΔt 
ηZ,k0 0dω2Jωω2sinωΔt  sin14ωΔt   coth12ωβ  cosk14ωΔt+isink14ωΔt 
ηZ,N,k 0dω2Jωω2sinωΔt  sin14ωΔt   coth12ωβ  cosNk14ωΔt+isinNk14ωΔt 
ηZ,k1 0dω2Jωω2sinωΔt  sin14ωΔt   coth12ωβ  cosk54ωΔt+isink54ωΔt 
ηZ,N−1,k 0dω2Jωω2sinωΔt  sin14ωΔt   coth12ωβ  cosNk54ωΔt+isinNk54ωΔt 
ηZ,kk 0dω2Jωω2sin2ωΔt     coth12ωβ  cos2kkωΔt+isin2kkωΔt 

Finally, we note that just as in the case of a single system coupled to a harmonic bath,60 the MPI influence functional coefficients could also be expressed in terms of correlation functions obtained from classical molecular dynamics calculations.

In this section, we illustrate the inclusion of harmonic degrees of freedom on a spin chain, consisting of two-level systems coupled to model bath described by a continuous spectral density, and a molecular aggregate, where we treat explicitly all available normal mode vibrations of each monomer.

Models of spin arrays exhibit intriguing dynamical effects, as well as quantum phase transitions. They are also of interest to magnetic materials, where external fields can be used to tune the spin tunneling parameter and thus alter the system’s dynamical behavior. Here, we explore the real-time dynamics of one-dimensional Heisenberg spin chains57 with nearest neighbor interactions (see Fig. 3) in the presence of dissipative interactions arising from coupled harmonic bath terms. The total Hamiltonian has the form

(5.1)

where σ^iα(i=x,y,z) are the Pauli spin matrices for unit α, Jy and Jz are the spin–spin coupling strengths, and E is the external field. The dissipative bath is characterized by an Ohmic spectral density,15 given by

(5.2)

where ωc is the cutoff frequency and ξ is the dimensionless Kondo parameter.

FIG. 3.

A one-dimensional chain of coupled spins. Each spin is coupled to a dissipative harmonic bath.

FIG. 3.

A one-dimensional chain of coupled spins. Each spin is coupled to a dissipative harmonic bath.

Close modal

Figure 4 shows the dynamics of a Heisenberg spin chain with r = 31 units, where each spin is coupled to its neighboring units, for two sets of spin–spin coupling parameters Jy=0.1E,Jz=0.5E and Jy=0.5E,Jz=0.5E. Even in the absence of coupling to a bath, calculating the dynamics of the Heisenberg chain presents a serious challenge. With the chosen value of r = 31 two-level systems, the electronic Hilbert space is equal to 231. Thus, the dissipative Heisenberg chain with 31 units is equivalent to a system of 231 (not just 31) states coupled to a bath. The bath parameters are set to ξ = 0.2, ωc=2E, and the temperature is kBT=E. We start with all spins aligned in the “up” state and show the probability Pα(t) for one of the edge spins (α = 1) to remain in the “up” state as a function of time, while the rest of the spins are traced over. The results were obtained using 2(N + 1) = 22 path integral variables for each nonterminal spin. The calculation for the last time point took about 77 min for each nonterminal linking step on a single AMD Epyc 7551P 2.0Ghz processor, with a total of approximately 36 h for the 31-unit chain.

FIG. 4.

Probability for the edge spin to remain in the “up” state in a Heisenberg chain of 31 spins for two sets of nearest neighbor couplings. Red lines with small markers show the dynamics of the chain when each unit is coupled to an Ohmic bath at kBT=E. Blue lines with small markers show the probability in the absence of dissipative bath degrees of freedom. Black solid lines show the dynamics of an isolated spin coupled to the same Ohmic bath. Black dashed lines show the tunneling oscillations of a single, isolated spin. Left panel: Jy=0.1E,Jz=0.5E. Right panel: Jy=0.5E,Jz=0.5E.

FIG. 4.

Probability for the edge spin to remain in the “up” state in a Heisenberg chain of 31 spins for two sets of nearest neighbor couplings. Red lines with small markers show the dynamics of the chain when each unit is coupled to an Ohmic bath at kBT=E. Blue lines with small markers show the probability in the absence of dissipative bath degrees of freedom. Black solid lines show the dynamics of an isolated spin coupled to the same Ohmic bath. Black dashed lines show the tunneling oscillations of a single, isolated spin. Left panel: Jy=0.1E,Jz=0.5E. Right panel: Jy=0.5E,Jz=0.5E.

Close modal

In the absence of dissipative interactions, an isolated spin undergoes coherent tunneling oscillations between the “up” and “down” states. As is well known,15 interaction of a single spin with a dissipative bath quenches the coherent dynamics, leading to underdamped oscillations or even monotonic decay to equilibrium. Figure 4 shows that spin–spin couplings are quite effective in quenching the amplitude of the single-spin oscillations, even in the absence of a dissipative harmonic bath. However, the resulting time evolution differs significantly from spin–boson dynamics. Furthermore, the two sets of spin–spin coupling parameters are seen to produce quite different dynamics. The first set of parameters with a smaller Jy coupling leads to an interesting interference effect, which causes an increase in the probability after the first recurrence of diminished amplitude. Interference is almost wiped out in the second case, where a blue shift in the oscillation of the probability is observed. When each spin is also coupled to a dissipative bath, the tunneling oscillations become further damped, and the interference characteristicsobserved with the first set of parameters disappears. Thus, we see that while coupling to neighboring spins generally leads to damping effects similar to those from harmonic dissipative baths, effects of spin–spin interactions are more complex and subtle, and the resulting dynamics depends delicately on the type and strength of coupling between spins.

In Fig. 5, we show the effects of increasing the temperature of the bath for the Heisenberg chain with Jy=0.1E,Jz=0.5E. As expected, with the increase in available thermal energy, the damping effects become more pronounced, and eventually, the oscillatory pattern almost disappears.

FIG. 5.

Survival probability of the “up” state of spin 1 in the Heisenberg chain with parameters from the left panel of Fig. 4, coupled to an Ohmic bath at different temperatures. Blue curve: kBT=0.2E. Violet curve: kBT=E. Red curve: kBT=5E.

FIG. 5.

Survival probability of the “up” state of spin 1 in the Heisenberg chain with parameters from the left panel of Fig. 4, coupled to an Ohmic bath at different temperatures. Blue curve: kBT=0.2E. Violet curve: kBT=E. Red curve: kBT=5E.

Close modal

As a second example, we apply the procedure described in Secs. II–IV to simulate the excitation energy transfer dynamics in J-aggregates of perylene bisimide (PBI; see Fig. 6), whose photochemical properties have been the subject of the recent electronic structure calculations.61,62 Each PBI molecule has Mα = 2 electronic states, nα = 1, 2, and its electronic Hamiltonian is described by Eq. (2.3) with hα,11e=hα,12e=hα,21e=0 and hα,22e=2.13 eV.62 Each electronic state is coupled to 28 intramolecular normal modes according to Eq. (2.4) with sα,1 = 0 and sα,2 = 2.The neighboring units interact via Frenkel exciton58 coupling terms,

(5.3)

The electronic energies, inter-monomer couplings, normal mode frequencies, and exciton-vibration coupling coefficients were obtained from electronic structure calculations.62 The spectral density is shown in Fig. 6 in terms of Huang–Rhys factors, which are related to the coupling coefficients through the expression

(5.4)
FIG. 6.

Huang–Rhys factors for the 28 intramolecular vibrational modes of perylene bisimide.62 For clarity, the lines for the two most strongly coupled modes have been reduced in height by a factor of 4.

FIG. 6.

Huang–Rhys factors for the 28 intramolecular vibrational modes of perylene bisimide.62 For clarity, the lines for the two most strongly coupled modes have been reduced in height by a factor of 4.

Close modal

Figure 7 shows the excited state population of the initially excited molecule in a small aggregate of r = 3 PBI molecules and a long aggregate containing r = 25 units at 300 K, with the initial excitation placed on the central unit. In the larger aggregate, a total of 700 discrete normal mode vibrations were treated explicitly via the fully quantum mechanical MPI procedure described in this paper. At the given temperature, the vibrational partition function of each molecular unit is equal to 2.8 × 1012, and thus, explicit treatment of all intramolecular modes via wavefunction-based methods is not practical. The calculation for the longest time point with N = 12 (26 beads) took about 10 min for each nonterminal linking step on a single AMD Epyc 7551P 2.0Ghz processor, with a total of approximately 220 min for the 25-unit aggregate. The coupling to molecular vibrations strongly affects the population dynamics, causing rapid damping of the electronic state population and some smearing of the minor recurrences observed in the long chain.

FIG. 7.

Survival probability of the initially excited state for PBI chains of (left) 3 and (right) 25 molecular units. Red curve: full electron-vibration dynamics. Black curve: dynamics in the absence of vibrations.

FIG. 7.

Survival probability of the initially excited state for PBI chains of (left) 3 and (right) 25 molecular units. Red curve: full electron-vibration dynamics. Black curve: dynamics in the absence of vibrations.

Close modal

The MPI representation of real-time dynamics enables the simulation of extended systems characterized by local interactions with effort that scales linearly with the number of units, allowing the study of long systems that are not amenable to wavefunction-based schemes. As a path integral method, MPI allows the treatment of any number of harmonic bath degrees of freedom that may couple to each unit, at zero or finite temperature, with practically no additional cost. This is possible by obtaining the influence functional for each unit in terms of pre-computed coefficients, which are specific to the particular discretization of the force on the bath arising from the particular splitting of the short-time propagators. In the common case where the coupling between adjacent units is diagonal in the given basis, the influence functional coefficients are identical to those in the QuAPI expression.25,63 When the coupling between units has a general, non-diagonal form, two-unit propagators must be employed with the swapped-pair ordering in order to avoid the unnecessary doubling of path integral variables. This ordering leads to different force sequences and thus different expressions for the discretized influence functional coefficients. In this paper, we have described the procedure for including the effects of harmonic bath modes to the MPI algorithm and derived the expressions required to obtain the influence functional coefficients for the general form of the MPI Hamiltonian. The method can be implemented with a model spectral density or with discrete bath modes representing molecular vibrations.

The methodology we have developed allows the simulation of time-dependent properties in extended systems representing a variety of different physical situations. For example, it can be used to obtain dynamical properties of spin arrays described by quantum Ising or Heisenberg Hamiltonians in which the spins are subject to dissipative effects. Even in the absence of dissipative interactions, these systems continue to present significant challenges to numerically exact methods, as the size of the Hilbert space grows exponentially with chain length. Thus, the MPI methodology can contribute to the understanding of many interesting physics of these systems, particularly the influence of dissipative phonon environments at finite temperature on the dynamics of spin arrays and their behavior relevant to quantum phase transitions. We emphasize that the MPI methodology is not restricted to simple linear arrays of units but can be extended to obtain the dynamics of more complex structures, such as branched chains, cyclic geometries, and dendrimers. We have recently presented MPI calculations of excitation energy transfer in the B850 ring of the light harvesting complex encountered in photosynthetic bacteria.64 

Another very interesting situation arises in the dynamics of molecular aggregates, which involve electronically coupled molecular units, each of which includes several tens or hundreds of intramolecular vibrational degrees of freedom. The resulting Hilbert space of these systems can be vast, particularly under ambient conditions, which lead to several populated vibrational states in many of the vibrational modes. The MPI methodology can treat such aggregates fully quantum mechanically. We have presented MPI calculations on excitation energy transfer in J-aggregates of perylene bisimide and (in a separate work64) investigated the energy transfer process in bacteriochlorophyll chains of varying length. Understanding such properties is important for designing structures for efficient solar energy harvest.

Finally, the MPI methodology can be used in conjunction with generalized quantum master equation65–67 (GQME) approaches in order to extend the time length that can be reached solely via the path integral calculation. GQME calculations require input obtained via different methods over the time span of the memory kernel.68–71 If this input is computed via numerically exact methods, then the GQME results at longer times will also be exact. Our recent work using QuAPI and QCPI calculations to build the GQME kernel,72 as well as the SMatPI decomposition35,36 of system–bath dynamics, suggests that the span of the GQME kernel is easily accessible by the time length of converged MPI calculations. Unlike the cases of QuAPI and QCPI, whose iterative decomposition is more efficient than GQME time-extensions,72 the combination of MPI with GQME is very appealing for obtaining numerically exact results on long molecular aggregates over long simulation times. Work along these lines will be reported in future papers by our group.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This material is based upon work supported by the Air Force Office of Scientific Research under AFOSR Award No. FA9550-18-1-0291. The calculations were performed on a computer cluster funded by the National Science Foundation under Award No. CHE-1665281.

1.
R. P.
Feynman
,
Rev. Mod. Phys.
20
,
367
387
(
1948
).
2.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
New York
,
1965
).
3.
N.
Metropolis
,
A. W.
Rosenbluth
,
M. N.
Rosenbluth
,
A. H.
Teller
, and
E.
Teller
,
J. Chem. Phys.
21
,
1087
1092
(
1953
).
4.
R. P.
Feynman
,
Statistical Mechanics
(
Addison-Wesley
,
Redwood City
,
1972
).
5.
E. L.
Pollock
and
D. M.
Ceperley
,
Phys. Rev. B
30
(
5
),
2555
2568
(
1984
).
6.
B.
Bernu
and
D. M.
Ceperley
, in
Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms
, edited by
J.
Grotendorst
,
D.
Marx
, and
A.
Muramatsu
(
John von Neumann Institute for Computing
,
2002
), Vol. 10.
7.
D. M.
Ceperley
,
Rev. Mod. Phys.
67
,
279
355
(
1995
).
8.
D.
Chandler
and
P. G.
Wolynes
,
J. Chem. Phys.
74
,
4078
4095
(
1981
).
9.
M.
Parrinello
and
A.
Rahman
,
J. Chem. Phys.
80
,
860
867
(
1984
).
10.
D.
Marx
and
M.
Parrinello
,
J. Chem. Phys.
104
,
4077
4082
(
1996
).
11.
M. E.
Tuckerman
,
D.
Marx
,
M. L.
Klein
, and
M.
Parrinello
,
J. Chem. Phys.
104
,
5579
5588
(
1996
).
12.
P. A.
Deymier
,
K.
Runge
,
K.-D.
Oh
, and
G. E.
Jabbour
, in
Multiscale Paradigms in Integrated Computational Materials Science and Engineering
, edited by
P. A.
Deymier
,
K.
Runge
, and
K.
Muralidharan
(
Springer
,
2015
), Vol. 226, pp.
13
106
.
13.
J. D.
Doll
,
D. L.
Freeman
, and
T. L.
Beck
,
Adv. Chem. Phys.
78
,
61
127
(
1990
).
14.
A. O.
Caldeira
and
A. J.
Leggett
,
Physica A
121
,
587
616
(
1983
).
15.
A. J.
Leggett
,
S.
Chakravarty
,
A. T.
Dorsey
,
M. P. A.
Fisher
,
A.
Garg
, and
W.
Zwerger
,
Rev. Mod. Phys.
59
,
1
85
(
1987
).
16.
17.
P.
Hänggi
,
P.
Talkner
, and
M.
Borkovec
,
Rev. Mod. Phys.
62
,
251
341
(
1990
).
18.
N.
Makri
,
J. Phys. Chem.
103
,
2823
2829
(
1999
).
19.
R. A.
Marcus
and
N.
Sutin
,
Biochim. Biophys. Acta
811
,
265
322
(
1985
).
20.
R.
Lambert
and
N.
Makri
,
J. Chem. Phys.
137
,
22A552
(
2012
).
21.
P. L.
Walters
and
N.
Makri
,
J. Phys. Chem. Lett.
6
,
4959
4965
(
2015
).
22.
W. H.
Miller
,
N. C.
Handy
, and
J. E.
Adams
,
J. Chem. Phys.
72
,
99
112
(
1980
).
23.
R. P.
Feynman
and
F. L.
Vernon
,
Ann. Phys.
24
,
118
173
(
1963
).
24.
N.
Makri
and
D. E.
Makarov
,
J. Chem. Phys.
102
,
4600
4610
(
1995
).
25.
N.
Makri
and
D. E.
Makarov
,
J. Chem. Phys.
102
,
4611
4618
(
1995
).
26.
J.
Shao
and
N.
Makri
,
Chem. Phys.
268
,
1
10
(
2001
).
27.
J.
Shao
and
N.
Makri
,
J. Chem. Phys.
116
,
507
514
(
2002
).
28.
M. M.
Sahrapour
and
N.
Makri
,
J. Chem. Phys.
132
,
134506
(
2010
).
29.
E.
Sim
and
N.
Makri
,
Chem. Phys. Lett.
249
,
224
230
(
1996
).
30.
E.
Sim
,
J. Chem. Phys.
115
,
4450
4456
(
2001
).
31.
M.
Richter
and
B. P.
Fingerhut
,
J. Chem. Phys.
146
,
214101
(
2017
).
32.
N.
Makri
,
Chem. Phys. Lett.
593
,
93
103
(
2014
).
33.
N.
Makri
,
J. Chem. Phys.
146
,
134101
(
2017
).
34.
A.
Strathearn
,
P.
Kirton
,
D.
Kilda
,
J.
Keeling
, and
B. W.
Lovett
,
Nat. Commun.
9
,
3322
(
2018
).
35.
N.
Makri
,
J. Chem. Phys.
152
,
041104
(
2020
).
36.
N.
Makri
,
J. Chem. Theory Comput.
16
,
4038
4049
(
2020
).
37.
H.-D.
Meyer
,
U.
Manthe
, and
L. S.
Cederbaum
,
Chem. Phys. Lett.
165
,
73
78
(
1990
).
38.
M. H.
Beck
,
A.
Jäckle
,
G. A.
Worth
, and
H.-D.
Meyer
,
Phys. Rep.
324
,
1
105
(
2000
).
39.
G. A.
Worth
,
H.-D.
Meyer
,
H.
Köppel
,
L. S.
Cederbaum
, and
I.
Burghardt
,
Int. Rev. Phys. Chem.
27
(
3
),
569
606
(
2008
).
40.
H.
Wang
and
M.
Thoss
,
J. Chem. Phys.
119
(
3
),
1289
1299
(
2003
).
41.
W.
Popp
,
M.
Polkehn
,
K. H.
Hughes
,
R.
Martinazzo
, and
I.
Burghardt
,
J. Chem. Phys.
150
(
24
),
244114
(
2019
).
42.
H.
Wang
,
J. Chem. Phys.
113
,
9948
9956
(
2000
).
43.
A.
Ishizaki
and
Y.
Tanimura
,
J. Phys. Soc. Jpn.
74
(
12
),
3131
3134
(
2005
).
44.
H.
Liu
,
L.
Zhu
,
S.
Bai
, and
Q.
Shi
,
J. Chem. Phys.
140
(
13
),
134106
(
2014
).
45.
46.
U.
Schollwöck
,
Ann. Phys.
326
,
96
192
(
2011
).
47.
R.
Orús
,
Ann. Phys.
349
,
117
158
(
2014
).
48.
J.
Ren
,
Z.
Shuai
, and
G.
Kin-Lic Chan
,
J. Chem. Theory Comput.
14
(
10
),
5027
5039
(
2018
).
49.
A. D.
Somoza
,
O.
Marty
,
J.
Lim
,
S. F.
Huelga
, and
M. B.
Plenio
,
Phys. Rev. Lett.
123
,
100502
(
2019
).
50.
N.
Makri
,
J. Chem. Phys.
148
,
101101
(
2018
).
51.
N.
Makri
,
J. Chem. Phys.
149
,
214108
(
2018
).
52.
S.
Kundu
and
N.
Makri
,
J. Chem. Phys.
151
,
074110
(
2019
).
53.
S.
Kundu
and
N.
Makri
,
Mol. Phys.
(submitted).
54.
S.
Suzuki
,
J.
Inoue
, and
B. K.
Chakrabarty
,
Quantum Ising Phases and Transitions in Transverse Ising Models
, 2nd ed. (
Springer-Verlag
,
Berlin-Heidelberg
,
2013
).
55.
J.
Dziarmaga
,
Phys. Rev. Lett.
95
,
245701
(
2005
).
56.
M.
Heyl
,
A.
Polkovnikov
, and
S.
Kehrein
,
Phys. Rev. Lett.
110
,
135704
(
2013
).
57.
W.
Heisenberg
,
Z. Phys.
49
,
619
636
(
1928
).
58.
59.
N.
Makri
,
Chem. Phys. Lett.
193
,
435
444
(
1992
).
60.
T. C.
Allen
,
P. L.
Walters
, and
N.
Makri
,
J. Chem. Theory Comput.
12
,
4169
4177
(
2016
).
61.
D.
Ambrosek
,
H.
Marciniak
,
S.
Lochbrunner
,
J.
Tatchen
,
X.-Q.
Li
,
F.
Würthner
, and
O.
Kühn
,
Phys. Chem. Chem. Phys.
13
(
39
),
17649
17657
(
2011
).
62.
D.
Ambrosek
,
A.
Köhn
,
J.
Schulze
, and
O.
Kühn
,
J. Phys. Chem. A
116
(
46
),
11451
11458
(
2012
).
63.
N.
Makri
,
J. Math. Phys.
36
,
2430
2456
(
1995
).
64.
S.
Kundu
and
N.
Makri
, “
Real-time path integral simulation of exciton-vibration dynamics in light harvesting bacteriochlorophyll aggregates
,”
Proc. Natl. Acad. Sci. U. S. A.
(submitted).
65.
S.
Nakajima
,
Prog. Theor. Phys.
20
,
948
(
1958
).
66.
R.
Zwanzig
,
J. Chem. Phys.
33
,
1338
1341
(
1960
).
67.
H.
Mori
,
Prog. Theor. Phys.
33
,
423
(
1965
).
68.
A.
Montoya-Castillo
and
D. R.
Reichman
,
J. Chem. Phys.
144
,
184104
(
2016
).
69.
A.
Kelly
,
A.
Montoya-Castillo
,
L.
Wang
, and
T. E.
Markland
,
J. Chem. Phys.
144
,
184105
(
2016
).
70.
L.
Kidon
,
H.
Wang
,
M.
Thoss
, and
E.
Rabani
,
J. Chem. Phys.
149
,
104105
(
2018
).
71.
E.
Mulvihill
,
X.
Gao
,
Y.
Liu
,
A.
Schubert
,
B. D.
Dunietz
, and
E.
Geva
,
J. Chem. Phys.
151
,
074103
(
2019
).
72.
S.
Chatterjee
and
N.
Makri
,
J. Phys. Chem. B
123
,
10470
(
2019
).