In a recent communication [N. Makri, J. Chem. Phys. **148**, 101101 (2018)], it was shown that the locality of interactions in many systems of interest allows a decomposition of the path integral and its evaluation via sequential linking of the paths of relatively small “modules” (e.g., chemical groups or monomers). The present paper describes the modular path integral methodology for simulating dynamical properties by propagating the density matrix in real time. The procedure is first presented for the simple topology of a single-file arrangement of units interacting via nearest neighbor couplings and subsequently extended to the calculation of two-particle correlations in arrays that may also contain some long-range interactions, to the treatment of systems with side chains or cyclic structures, to the simulation of internal dynamics in long organic molecules, and to the modifications required for coupling of one or several units of a system to dissipative environments. Illustrative applications to the dynamics of interacting two-level-systems are presented.

## I. INTRODUCTION

Quantum mechanical calculations in polyatomic molecules or condensed phase systems are plagued by the exponential growth of wavefunction storage with system size. Overcoming the exponential scaling of quantum mechanics without resorting to approximations has been possible in some situations, and generally exploits one of two avenues: stochastic sampling or tensor decompositions. In connection with the path integral^{1,2} representation of quantum statistical mechanics,^{3} Monte Carlo methods^{4} offer a powerful, numerically exact treatment to the equilibrium properties of many-particle (spinless or boson) systems.^{5–7} The quantum-classical isomorphism of the imaginary time path integral^{8} also allows the use of molecular dynamics algorithms as an alternative to Monte Carlo sampling.^{9} These stochastic methods are excellent at estimating equilibrium properties, as the integrand arising from the Boltzmann operator is smooth and localized in space. Unfortunately, extension of these methods to real-time properties has not met with success because of the “sign problem,” i.e., the dramatic phase cancellation associated with the oscillatory real-time propagator.^{10,11}

Tensor decompositions of the wavefunction have emerged as accurate and efficient tools for studying some types of strongly correlated systems. In particular, the density matrix renormalization group^{12} (DMRG) approach has led to powerful algorithms for computing ground state properties of extended one-dimensional systems, although finite-temperature^{13} and time-dependent^{14} variants are also available. Furthermore, tensor decompositions *in time* of the real-time path integral^{15–17} have led to the formulation of numerically exact algorithms for simulating the dynamics of quantum dissipative systems that scale linearly with propagation time.

A recent communication^{18} introduced a tensor decomposition of the path integral *in space*, which appears promising for computing dynamical properties in large systems characterized by some degree of locality, i.e., where potential interactions of each particle or a small group of particles (a “module”) extend mostly to neighboring modules. Instead of propagating in time a spatial representation of the many-particle wavefunction or density matrix, the modular path integral (MPI) algorithm proceeds by sequentially linking the paths of adjacent modules and leads to linear or (since identical modules are treated just once) sublinear scaling with system size. The basic requirement of the MPI approach is the ability to store the Feynman paths. For processes that reach equilibrium in a relatively short time (e.g., 10-20 path integral time steps), storage of the paths is certainly within reach. Longer time calculations can utilize stochastic sampling and weight-based filtering ideas^{19–23} employed in the context of quantum dissipative systems. The first illustration of the MPI decomposition^{18} focused on complex-time correlation functions in arrays of 50 coupled spins.

The present paper describes the MPI methodology for the real-time propagation of density matrices in a variety of situations. Section II presents the algorithm for a simple, linear arrangement of units connected via nearest-neighbor interactions. A number of extensions are described in Sec. III. These include systems of different topologies such as branched chains and cyclic structures, the use of atom-specific path integral time steps, the possibility of including some long-range interactions, the calculation of two-particle correlations, the simulation of intramolecular dynamics in long organic molecules, and the inclusion of dissipative effects. Model applications to arrays of two-level systems (TLSs) (e.g., spins or coupled excitons) are given in Sec. IV, and some concluding remarks are given in Sec. V.

## II. THE BASIC ALGORITHM: DENSITY MATRIX FOR A SIMPLE ARRAY OF MULTILEVEL SYSTEMS

Consider a system composed of units (or monomers) labeled *α* = 1, 2, …, *n*, each of which is defined in terms of the *M*_{α} discrete states $\varphi i(\alpha )$, such that the total Hamiltonian has the form

where

are one-particle operators and the term

couples each unit to its nearest neighbors. The discrete states might represent “sites,” as in the case of the left- and right-localized states of a two-level system (TLS), or they may have been obtained from a discrete variable representation^{24} (DVR) of a continuous Hamiltonian. In either case, it is assumed that the potential-like operators *V*^{(α)} are local in the basis of these discrete states, and that they couple the units to their nearest neighbors. Such discrete states form the optimal representation for grid-based^{25} real-time path integral methods,^{26} eliminating the rapid oscillations of the real-time propagator^{27} and offering the possibility of quadrature evaluation of the path sum. Note that the number *M*_{α} of states in each unit, as well as the intra-unit Hamiltonian and monomer-to-monomer coupling parameters, may be different for each unit.

The form of the Hamiltonian in Eq. (2.1) is relevant to a plethora of physical situations. For example, it describes Heisenberg and quantum Ising chains,^{28} which are important in the understanding of ferromagnetic materials and in the study of quantum phase transitions. With minor modifications, this Hamiltonian can be adapted to describe the vibrational dynamics of simple molecular chains (e.g., cumulenes), as well as the dynamics of exciton transfer in single file-arranged molecular aggregates.

Dynamical observables, at zero or finite temperature, can be obtained from the evolution of the density operator

where $Tr\u2009\rho ^(0)=1$. The initial condition may correspond to a pure state, or may describe a mixed ensemble, e.g., $\rho ^(0)$ may be given by the Boltzmann operator. Throughout this section, it is assumed that the initial density operator factorizes, i.e.,

Upon propagating the density operator in time, multiplying with an operator $O$ and taking the trace with respect to all units, one obtains the expectation value $O(t)$ of this operator as a function of time. The operator $O$ may correspond to a property of a particular unit, or to a property of a number of units. In this section, it is assumed that $O^=\u220f\alpha =1nO^(\alpha )$. For example, if a property of unit *α*_{1} is sought, $O^(\alpha )=\xce(\alpha )$ (the identity operator for monomer *α*) for all *α* ≠ *α*_{1}.

The path integral representation of the density operator is obtained by expressing the time evolution operator (and its inverse) in terms of shorter time operators using the exponential identity

where Δ*t* = *t*/*N* is the path integral time step, using the resolution of identity

and applying the Trotter factorization,^{29}

With the Hamiltonian of Eq. (2.1), the Trotter-factored short time propagator becomes

The single-unit factors can be obtained by using the eigenstates and eigenvalues of the corresponding operators,^{25,27} i.e.,

where

Using Eq. (2.9), one arrives at the discretized path integral representation for the expectation value of $O$,

where $wj(n+1)\u22610$. For each unit, the sequence of states $\varphi i\alpha ,0+(\alpha ),\u2026,\varphi i\alpha ,N+(\alpha )$ forms a discrete path in the forward time direction, while the sequence $\varphi i\alpha ,0\u2212(\alpha ),\u2026,\varphi i\alpha ,N\u2212(\alpha )$ defines a path along the backward time direction. Thus, there are (*M*_{α})^{2N+2} forward-backward paths for each monomer, and Eq. (2.12) contains a total of $\u220f\alpha =1nM\alpha 2N+2$ terms. The instantaneous states along the paths are the discrete analogue of the path coordinates or “beads,” which are encountered in equilibrium path integral calculations in continuous coordinate space. The use of discrete states in the path integral (made possible via a DVR of the path integral,^{26} even in the presence of an influence functional) has been shown to be the optimal representation for real-time calculations.

Even with the smallest possible value of *M*_{α} = 2, the number of terms is 2^{(2N+2)n}, which typically is not amenable to direct summation. As discussed in the Introduction, Monte Carlo methods encounter a disastrous cancellation in real-time path integral calculations.^{10,11,30} The most promising approach is to look for decompositions of the integrand, which would allow evaluation of the sum in smaller (thus manageable) tasks involving ordinary quadrature or a combination of quadrature and Monte Carlo procedures. The latter has proven to be a very successful strategy in the context of the quantum-classical path integral (QCPI),^{31–33} which has been shown feasible even in the presence of oscillatory phase factors from thousands of atoms.^{34}

To this end, examination of Eq. (2.12) reveals two types of interactions among the path integral variables. First, the path integral variable of unit *α* at the time *k*Δ*t* (for 0 < *k* < *N*) is coupled to the path integral variables of the same unit at the time points (*k* − 1)Δ*t* and (*k* + 1)Δ*t*. Second, this path integral variable is also linked to the path integral variables of the two adjacent atoms at the same time point *k*Δ*t*. This connectivity is illustrated in Fig. 1.

The conventional approach for wavefunction propagation of small systems, which is also applicable to the density matrix (with higher computational cost), involves storing the initial state (in the chosen grid or basis representation) andpropagating in time via matrix-vector multiplication.^{35} In the schematic illustration of Fig. 1, the initial condition is denoted by shaded rectangles and propagation would occur in the left and right directions. Even though the density matrices of the monomers are initially uncoupled (because of the current assumption of a factorizable density), they would become entangled after the first step of time propagation, and subsequent propagation would require storage on the set of states spanning the full Hilbert space of the multi-monomer system. This Hilbert space is prohibitively large in the case of a long chain (*n* ≫ 10).

Instead, the MPI scheme proceeds by forming all paths of the first unit and proceeding to add their amplitudes after linking them to those of the neighboring unit, i.e., performing the calculation in the vertical direction. While the storage of all single-unit paths is not always practical, there are many situations where this task is feasible. Many properties of interest are obtained from propagation to only *N* ≃ 10–20 time steps, in which case (at least for monomers consisting of 2-3 states) storage of the Feynman paths is relatively straightforward. Furthermore, there are several systematic methods for eliminating paths that make exponentially small contributions.^{19–23} Since the number of retained paths can amount to an exponentially small fraction of those possible, the huge gain in efficiency implies that propagation to considerably longer times should be possible. Additional enhancements and new possibilities for extending the MPI methodology to longer times will likely be identified.

The MPI procedure begins with the construction of all forward-backward paths of each unit (with initial and final states that conform to the structure of the initial density matrix and the observable operator), along with their amplitudes, and the matrices that link adjacent units in the path integral expression. For notational simplicity, each forward-backward path pair $\varphi i\alpha ,0+(\alpha ),\u2026,\varphi i\alpha ,N+(\alpha );\varphi i\alpha ,0\u2212(\alpha ),\u2026,\varphi i\alpha ,N\u2212(\alpha )$ is labeled by the unit number *α* and the sequence of state indices $i\alpha ,0\xb1,i\alpha ,1\xb1,\u2026,i\alpha ,N\xb1$. The single-monomer path amplitude is given by the product

The monomer linking matrices **T**^{(α+1,α)} couple each path of a monomer to a path of its nearest neighbor and are given by

The linking matrices can be computed as needed and do not need to be stored.

The MPI procedure is initialized by generating and temporarily storing the array **R**^{(1)} of path amplitudes for monomer 1,

Next, monomer 1 is linked to monomer 2 through multiplication by the linking matrix **T**^{(2,1)}. The paths of monomer 1 do not appear elsewhere in the path integral expression, so they can be summed as soon as they are linked to the paths of monomer 2. Once this summation has been performed, the path amplitudes of monomer 2 are included. This produces the array **R**^{(2)} of path amplitudes for monomer 2, augmented by its interaction with monomer 1. This sequence of steps is repeated to link monomer 2 to unit 3, and continued until the end of the chain is reached. The unit linking procedures correspond to the following expression, which is repeated to link each unit with the next unit on the chain,

where *α* = 1, …, *n* − 1. Equation (2.16) involves multiplication of the vector **R**^{(α)} by the square matrix **T**^{(α+1,α)}, followed by multiplication of each element of the resulting vector by the corresponding element of the vector **A**^{(α+1)}. Finally, once the end of the chain is reached, the observable of interest is obtained by summing the amplitudes with respect to all paths,

## III. MODIFICATIONS AND EXTENSIONS OF THE BASIC ALGORITHM

Several extensions of the procedure are possible and necessary to simulate a variety of systems. The most important of these are outlined here.

### A. Two-particle operators

Suppose one is interested in the correlations of monomers *α*_{1} and *α*_{2} at the time *N*Δ*t*, measured through an operator $Q(\alpha 1,\alpha 2)$, such that

To obtain the expectation value of $Q$, one must refrain from summing with respect to all states corresponding to the last time point of unit *α*_{1}, instead augmenting the propagated array with the coordinates of that time point,

Subsequent linking operations proceed as usual, until unit *α*_{2} is reached. At that point, the two-particle operator is inserted and the auxiliary coordinate of unit *α*_{1} is summed,

Since only a single additional coordinate is stored, the increase in required storage is minimal. The linked coordinates are illustrated by triangles in Fig. 1.

### B. Long-range interactions

Even though the efficiency of the MPI algorithm stems from the locality of couplings among units, it is possible to include some long-range interactions at increased computational cost. This is particularly useful for long molecules that (in addition to nonpolar carbons and hydrogens, which are usually described in terms of bond stretching, bending, and torsional interactions) include a pair of polar atoms which interact via Coulombic terms.

To account for interactions between the distant units *α*_{1} and *α*_{2}, it is necessary to extend the procedure described in Subsection III A to couple the paths of these units at all times. Suppose this interaction adds the term $V^(\alpha 1)V^(\alpha 2)$ to the Hamiltonian. Upon reaching the first of these units, the propagated array begins to be augmented by the coordinates of its paths

etc. Once monomer *α*_{2} is reached, its paths are linked to those of unit *α*_{1} through multiplication with the phase factors corresponding to the $V^(\alpha 1)V^(\alpha 2)$ coupling, in a way similar to the inclusion of the nearest neighbor interactions in the path linking matrix. This involves the new linking matrix between these non-nearest neighbor monomers,

The amplitudes are then summed over the paths of unit *α*_{1} and the auxiliary storage of these paths is discontinued,

This procedure is illustrated in the bottom right panel of Fig. 1.

### C. Use of unit-dependent time steps

Heavier atoms (or atoms connected via lower frequency potentials) may be treated via larger path integral time steps, provided the time steps of all particles are multiples of a common value. This is routinely done in path integral Monte Carlo (PIMC) calculations and can result in a significant increase of efficiency. In the case of MPI, the reduction in the number of beads for certain units can lead to even higher savings, as these units may contain more states that would give rise to a larger number of paths.

The modification of the MPI algorithm to account for time step doubling is entirely straightforward and analogous to the corresponding PIMC procedure. This modification is illustrated in the middle right panel of Fig. 1. If desired, the time step can be increased by other integer factors in a similar way.

### D. Coupling to dissipative heat baths

In many situations of interest, the chain may be attached to a solid on one or on both ends, whose phonons act as a dissipative environment. The interaction of unit *α* with the phonons of the solid is described quantitatively by coupling that unit to a harmonic dissipative bath,^{36} i.e., augmenting its Hamiltonian by the term

where

The path integral representation of the system-bath Hamiltonians leads to an influence functional^{37} that augments the amplitudes of the system. In its discretized form, the influence functional that augments the path amplitude of unit *α* is given by

where the coefficients $\eta k\u2032k\u2033(\alpha )$ are available^{16} in terms of frequency integrals that involve the spectral density function^{36} characterizing the particular bath. The influence functional coefficients can also be obtained from the force autocorrelation function of the bath.^{38} This is convenient for mapping a complex environment on a harmonic bath from information available through classical molecular dynamics methods. Inclusion of the influence functional involves simple multiplication of the amplitudes of each unit, i.e.,

Since the MPI scheme involves summing with respect to the entire paths, no memory truncation is required. The influence functional interactions are illustrated in the top panel of Fig. 1.

### E. Boltzmann density

In many situations, the initial density matrix corresponds to the canonical ensemble for the composite multi-unit chain, which is not factorizable because of intermonomer couplings. One can follow either one of the following procedures to adapt the MPI scheme to this situation.

The most straightforward way involves using the imaginary-time path integral representation of the Boltzmann operator to bring the chain to the desired temperature. The monomer-path connections are similar to those described in Sec. II, but the first segment of each path would now correspond to imaginary time. After the necessary number of short imaginary-time steps have been taken, the paths would continue in the real time direction. Observables are collected from the real-time portion of the paths. The procedure requires longer paths, which now include imaginary-time Boltzmann segments. This approach can be used to calculate equilibrium correlation functions.

Another alternative is motivated by the use of Langevin thermostats in classical molecular dynamics simulations.^{39} One starts with a simple, factorized initial density, and uses the MPI methodology to propagate the chain in contact with a dissipative bath in real time (without including measurement operators in the propagator). Because systems evolve to equilibrium, this process will generate the density at the given temperature.

### F. Branched chains and cyclic structures

To discuss how to deal with chain branching, it is easier to think in terms of two chains, labeled A and B, with *n*_{A} and *n*_{B} monomers, respectively, connected to the first unit of chain C. The connections are illustrated in Fig. 2.

Starting from unit 1 of chain A, one proceeds toward the direction of the joint, storing the array $R(A,nA)$ when its last unit is reached. Similarly, starting from monomer 1 of chain B, one proceeds toward the joint and generates the array $R(B,nB)$. The amplitudes of these two chains are linked to the first unit of chain C through the operation,

This operation involves two separate summations, which are performed sequentially. The increase in cost is linear with respect to the number of different side chains. In fact, once the internal dynamics of a group comprising several units has been dealt with (by linking the monomers), the group acts as a new, higher-level module, and calculating the dynamics of the composite system consists of connecting such multi-unit modules. Note also that identical side chains can be dealt with just once and stored. This is particularly useful in calculations where there are repeating units, e.g., functional groups in a long molecule, or identical chromophores in a molecular aggregate.

The MPI treatment of a cyclic structure is almost identical to the treatment of long-range interactions described in Sec. III B. The MPI array is augmented with the path amplitudes of unit 1, until the last unit is reached. At this point, the paths of the last unit are linked to those of unit 1 and the array is summed with respect to all paths.

### G. Organic molecules

Finally, consider the internal dynamics of an organic molecule described by a force field with bonded interactions. (Extension to include distant electrostatic interactions would proceed as described in Sec. III B.) As an example, Fig. 3 presents the first few atoms of a conjugated alkene. All atoms are numbered.

For clarity of presentation, suppose first that only bond stretching interactions are included in the potential. Stretching interactions involve two atoms, thus the basic module consists of a single atom. In this case, the molecule is topologically identical to the case of a branched chain discussed in Sec. III F, with each hydrogen atom comprising a side chain. The MPI algorithm would couple (separately) the path amplitudes of H_{1} and H_{2} on C_{1}, then sum over those paths. Next, the paths of H_{3} would be coupled to those of C_{2} and summed. The paths of C_{1} (augmented by the interactions with the first two hydrogen atoms, which have already been accounted for) would then be connected to those of C_{2} and summed. Continuing this way, one proceeds along the molecular backbone with effort that grows linearly with molecular size.

Next, suppose that bond bending interactions are included as well. Bond angles involve three atoms, so the basic module must be augmented such that directly interacting, adjacent-bond atoms belong to neighboring modules. This implies that the basic module must include two adjacent carbon atoms, along with the attached hydrogens. In the example of Fig. 3, this can be achieved by defining module 1 to include C_{1}, H_{1}, H_{2}, C_{2}, and H_{3}, module 2 to include C_{3}, H_{4}, C_{4}, and H_{5}, etc. While the MPI decomposition again leads to linear scaling, the required effort is now greater because of the larger number of atoms in the basic module. However, it is well known that high-frequency CH stretching vibrations do not couple appreciably to the remaining internal motions; thus, CH bond lengths are often constrained to remain fixed in molecular dynamics simulations. The quantum mechanical equivalent of constraining the CH lengths is to use just the ground state for the CH vibrations, thus maintaining a relatively small number of states in each module. Finally, allowing for torsional interactions is possible by extending these ideas.

## IV. APPLICATION

The MPI methodology is illustrated below with a model chain or two-level systems (TLS), i.e., *M*_{α} = 2, connected through Ising-type (nearest neighbor) interactions. The Hamiltonian is given by the expression

where $\sigma ^x$, $\sigma ^z$ are the Pauli spin operators and $2\u210f\Omega \alpha $ are the TLS tunneling splittings. This Hamiltonian has the form of Eqs. (2.1)–(2.3) with

where Ω_{α} are the TLS tunneling splittings and *c*_{α} is the strength of the coupling between adjacent TLSs. The TLS chain is a zero temperature and initially all TLSs are in site 1, i.e.,

The population of site 1 of various TLSs is followed as a function of time under various conditions.

To illustrate the effects of an external dissipative environment, the results are also presented with a single TLS attached to a harmonic bath. In this case, the MPI with dissipation is identical to the quasiadiabatic propagator path integral (QuAPI) methodology.^{15–17} The coupling is given by Eq. (3.8) where the TLS “position” coordinates are $s1(\alpha )=1$, $s2(\alpha )=\u22121$. The chosen bath is characterized by the common Ohmic spectral density,^{40}

where *ω*_{c} is the cutoff frequency and $\xi $ is the dimensionless Kondo parameter.^{36}

Figure 4 shows the site population of the first (*α* = 1), third, (*α* = 3), and an inner (*α* = 10) TLS in a chain of *n* = 40 TLSs with all *c*_{α} = Ω_{α} ≡ Ω. The MPI calculations were performed on a laptop and utilized up to 2*N* = 16 beads. It is seen that the coupling to the other TLSs of the chain induces dissipative behavior over the time length examined, and the TLS populations exhibit an apparent decay. (Of course, it is clear that any finite system will eventually exhibit recurrences, but with 40 degrees of freedom these are not observed until much longer times.) The edge TLSs, which are linked to the rest of the chain only on one side, maintain some coherence with these parameters and exhibit quenched oscillatory dynamics. Coherence is largely quenched for the inner TLSs, and the third TLS already displays incoherent decay.

Also shown in Fig. 4 are results for the same chain, but with the first TLS interacting with a dissipative harmonic bath characterized by *ξ* = 0.5 and *ω*_{c} = 2 at two temperatures, corresponding to $\u210f\Omega \alpha \beta =5$ and $\u210f\Omega \alpha \beta =1$, as well as the population dynamics of a single TLS interacting with this dissipative bath under similar conditions. In line with well-established behaviors,^{40,41} the single TLS exhibits weak oscillations that are quenched faster in the case of the higher-temperature environment. It is seen that (with the chosen parameters) even in the absence of thermal fluctuations, the dissipative effects of the 40-TLS chain on the first TLS are stronger than those of the finite-temperature harmonic dissipative baths, producing more incoherent dynamics. The additional dissipative effects of the harmonic bath on the TLS to which the bath is attached are substantial, with the higher temperature bath leading to faster quenching of the tunneling oscillation. The harmonic bath causes a slightly faster decoherence of the third TLS. The effects of the bath on more distant TLSs are negligible with the given parameters.

## V. CONCLUDING REMARKS

The MPI decomposition of the path integral offers a numerically exact, fully quantum mechanical approach to the dynamics of extended systems with mostly local interactions. As a path integral-based method, it is applicable to zero- or finite-temperature properties and may be used to calculate the time evolution of expectation values and time correlation functions. The MPI algorithm does not rely on periodicity or require a simple composition of the system under study, thus it is easily applicable to complex molecules with branched backbones and rings, as long as the potential interactions are relatively local. Furthermore, some distant interactions may be included at somewhat higher cost, and dissipative effects arising from interaction with phonons are straightforward to include.

Local interactions are prevalent in chemical and physical systems of interest, and the MPI approach should provide useful results for those systems which may not be obtainable by other methods. For example, one-dimensional arrays of interacting spins, which serve as important models for quantum phase transitions, are ideally suited to this approach. The specifics of intramolecular energy flow through long organic molecules, which can be adequately described in terms of spectroscopically determined force fields with local bonded-atom potential interactions, remain poorly understood, as anharmonic terms are important, while classical trajectory studies are unable to properly account for zero-point energy. By explicitly treating potential anharmonicity, while fully accounting for quantum mechanical effects, MPI calculations can shed light on such processes. Exciton transport in single-file arrangements of molecular aggregates, for instance those encountered in photosynthetic antenna complexes, should also be amenable to the MPI methodology. In addition, the MPI scheme can easily incorporate external driving interactions.

Finally, it should be possible to combine the MPI treatment of an excitonic array with the QCPI^{31–33} treatment of molecular vibrations, leading to a rigorous and highly accurate treatment of vibronic effects in the dynamics of these systems.

## ACKNOWLEDGMENTS

This material is based upon work supported by the Air Force Office of Scientific Research under AFOSR Award No. FA9550-18-1-0291.