Modular path integral for discrete systems with non-diagonal couplings

Many systems with interesting physical, chemical, or biological properties are characterized by a one-dimensional or quasi-one-dimensional topology and mostly short-range interactions. Typical examples are long organic molecules such as linear or branched hydrocarbons, some of which possess conducting properties, molecular aggregates with important light harvesting functions, and structures with ferromagnetic properties of interest in materials research. In spite of their quasi-one-dimensional topology, these systems are often characterized by a considerable amount of structural complexity and many interacting electronic states which are coupled to local molecular vibrations and/or phonons. This complexity leads to intriguing many-body equilibrium and dynamical effects, whose understanding hinges on the availability of accurate quantum mechanical simulation methods.

The infamous exponential scaling of full-space wavefunction-based methods typically restricts attention to relatively small systems. The density matrix renormalization group¹ (DMRG) method, which takes advantage of tensor factorization along with singular value decomposition, has led to accurate and efficient tools for studying the ground-state properties of extended one-dimensional systems, and finite-temperature² and time-dependent³ variants of DMRG are also available. However, such methods cannot account simultaneously for the important effects that result from the coupling of hundreds of molecular vibrations and phonons to the large number of relevant electronic states of conducting polymers or magnetic materials.

Feynman’s path integral formulation^4,5 of quantum statistical mechanics,⁶ combined with stochastic Monte Carlo sampling methods,⁷ offers a powerful, numerically exact approach to thermal equilibrium properties of many-particle spinless or bosonic systems described by continuous coordinates. Furthermore, the quantum-classical isomorphism of the imaginary time path integral⁸ also allows the use of molecular dynamics algorithms as an alternative to Monte Carlo sampling.⁹ However, application of these tools to systems containing identical fermions, or to systems characterized by some discrete Hamiltonians, often encounters a severe “sign problem” that prevents convergence. Perhaps even more importantly, path integral Monte Carlo methods are impractical for computing dynamical properties, as the oscillatory nature of the real-time propagator leads to a massive sign problem even in the case of spinless systems with continuous coordinates.^10,11 Until very recently, numerically exact non-Monte Carlo-based real-time path integral methods have been available only for system-bath Hamiltonians.^12,13

Unlike methods that are intimately tied to wavefunctions, the path integral approach offers the possibility of treating normal mode vibrations and phonons fully quantum mechanically, simultaneously with the treatment of the coupled electronic states. Our group has invested much effort in this direction, mainly through the development of the iterative quasi-adiabatic propagator path integral (i-QuAPI) methodology^14–18 for evaluating the path integral with the Feynman-Vernon influence functional¹⁹ and its more recent extension that accelerates convergence in the most challenging strongly incoherent regime.^20,21 Recent work has extended these capabilities to processes in anharmonic environments (such as those from a solvent or protein, treated at the all-atom level) through the development of the quantum-classical path integral^22–27 (QCPI), a rigorous quantum-classical formulation that allows simulations of nonadiabatic dynamics without ad hoc assumptions or adjustable parameters.²⁸ These methods offer the possibility of simulating a number of important processes (such as electron transfer, barrier crossing, etc.) in a condensed phase or biological environments with unprecedented accuracy and are not limited by assumptions on the type of potential interactions among the nuclei comprising the system’s environment. Yet, the numerical effort for evaluating the path integral grows exponentially with the number of system states (or sites), and while various filtering schemes^20,29–31 can dramatically reduce this cost, path integral simulations of systems with hundreds of states remain out of reach.

The present paper continues the recent development of the modular path integral^32,33 (MPI) decomposition, which shows promise for treating both equilibrium and dynamical effects in systems with a quasi-one-dimensional topology and mostly local interactions, such as those discussed in the first paragraph. Rather than propagating in time a spatial representation of the many-particle wavefunction or density matrix, the MPI algorithm proceeds by sequentially linking the Feynman paths of adjacent units (the “modules”) and leads to linear (or, since identical modules are treated just once, sublinear) scaling with the system size, regardless of the exponential growth of the full Hilbert space.

The original MPI papers presented the algorithm for simple and branched chains of discrete systems with interactions between pairs of units that involved operators which are diagonal in a product basis of single-unit states. Obtaining a path integral representation of the propagator was straightforward in that case, and the main goal of these papers was to develop an efficient MPI decomposition. The present paper extends the MPI formulation to Hamiltonians where no product basis of single-particle states diagonalizes the coupling operator. This situation is often encountered in spin models, as well as the Frenkel exciton Hamiltonian.³⁴ The problem associated with nondiagonal operators is finding a convenient set of path integral variables, which allows an efficient MPI decomposition.

In Sec. II, we consider various ways of constructing propagators for Hamiltonians with intermonomer couplings of a general form and discuss their various features and drawbacks in the context of the MPI decomposition. By examining the diagrammatic representations of these propagators, we are able to identify in Sec. III a form that leads to a path integral representation that satisfies the detailed balance property while minimizing the number of variables and allowing optimal MPI scaling. In Sec. IV, we describe the MPI decomposition and its factorization. We illustrate the method in Sec. V, with application to two systems of coupled two-level systems (TLSs), which are motivated from interacting spin models and exciton energy transfer. Some concluding remarks, along with the forthcoming extension of the MPI methodology to account for molecular vibrations, are given in Sec. VI.

We consider an array of r monomers, described by a Hamiltonian of the form

where the superscript α = A, B, ... is the monomer index and $ĥα$ are single-monomer Hamiltonians with M_α discrete states, and the coupling between adjacent monomers is given by Hermitian operators that have the form

Here, $nα$ are basis states of monomer α, which may or may not be eigenstates of the single-monomer Hamiltonians.

For convenience, we label the individual monomers sequentially in the alphabetical order, so that the entire system has the structure A-B-C-D-…. We now construct pairs of adjacent units by dividing the Hamiltonian,

where (assuming there are more than four monomers)

etc.

Common Hamiltonians in continuous (e.g., Cartesian or normal mode) coordinates involve a potential coupling that is a function of position coordinates (which may be replaced by a discrete variable representation³⁵). Because the potential operator is local in the position space, the simple Trotter splitting³⁶ of the time evolution operator into single-particle and coupling factors allows the use of position states as auxiliary path integral variables, with the single-particle Hamiltonians entering in the form of short-time propagators that are easily evaluated numerically without further approximations.¹⁴ While imaginary time path integral calculations are typically performed in continuous coordinates and path integral variables are sampled by stochastic methods, the need to avoid the Monte Carlo sign problem in the case of real time evolution necessitates the use of discrete variables. System-dependent discrete variable representations (DVR) of the path integral³⁷ allow the use of such propagator factorizations within minimal-sized basis sets that may be chosen from system eigenstates³⁸ or a moving basis.³⁹ These discrete representations of the path integral should allow a straightforward MPI decomposition of systems with continuous potentials. The simple Trotter-like factorization of interparticle coupling is also possible for some discrete models, e.g., in the quantum Ising model, where the spin-spin interaction is a product of single-spin Pauli matrices and thus is local in a product basis of single-particle states. In all these cases, the interparticle coupling of Eq. (2.2) involves products of single-unit components and thus is diagonalizable in a monomer product basis.

However, there are many situations with two- or multistate units, where the coupling operator in Eq. (2.2) is not diagonal in any basis that is a product of single-monomer states. Famous examples in this category are the Frenkel exciton Hamiltonian,³⁴ which is used to describe energy transfer in conjugated polymers and chromophore arrays, and various models of interacting spins, which are important in ferromagnetic materials. In these situations, it is not possible to find path integral variables in which the short-time propagator is a product of single-particle components.

In the absence of a single-particle product basis that diagonalizes the coupling, we express the path integral in the product basis of the single-monomer states $nα$ and seek efficient ways of obtaining short-time propagators for use in the path integral. We note that the choice of these single-particle basis states is arbitrary. We choose basis states motivated by the physical process of interest and the relevant boundary conditions. For example, in the case of Frenkel exciton transfer, the basis states are the ground and excited eigenstates of the single-monomer Hamiltonian. When dealing with spin models, it is natural to use the site representation, i.e., the localized “right/left” or “up/down” states as the basis. The advantage of choosing the basis states in this fashion is that the boundary conditions of interest involve a single state (e.g., all spins start in the “up” configuration or a single monomer is excited initially with all other monomers in the ground state).

The first step is to obtain a short-time propagator for the Hamiltonian of Eq. (2.1) in a product basis, i.e., we seek to evaluate the propagator elements,

where Δt is the path integral time step, and the subscripts 1 and 2 correspond to the time points t₁ and t₂ = t₁ + Δt. We proceed to factor the time evolution operator as follows:

where $ÛABΔt=e−iĤABΔt/ℏ$⁠, etc. Equation (2.6) is an asymmetric Trotter factorization, where the leading error is proportional to the sum of nonvanishing commutators, $[ĤAB,ĤBC]+⋯Δt2$⁠, etc. While there are many orderings possible in Eq. (2.6), we note that reversing the order of the factors,

produces a Trotter factorization where all commutators have the opposite sign. (It is obvious that $ÛAB=ÛBA$⁠.) Thus, by using the arithmetic mean of Eqs. (2.6) and (2.7), one obtains a short time propagator where the quadratic term in Δt vanishes, producing a symmetric factorization with an error of order Δt³. (We note parenthetically that in the case of continuous potentials, the scaling of the Trotter factorization error with the time step, which is rigorously derived via expansion of the time evolution operator, does not apply to the propagator in the continuous coordinate representation,⁴⁰ even though symmetrization allows larger time steps even in that case.)

To obtain the propagator in the product basis of single monomer eigenstates, it is necessary to insert the resolution of identity in terms of auxiliary variables between adjacent monomers,

The first and last monomers do not require auxiliary variables. Each auxiliary variable appears together with a state indexed 1 in the kets, but the same variable is placed next to a state indexed 2 in the bras. Thus, the auxiliary variables do not correspond to the time associated with either of the two propagator endpoints, so we label them with the subscript 12. The reverse-order factorization of the short time evolution operator leads to the propagator expression

The connections among the propagator endpoint and auxiliary variables are illustrated in Fig. 1. Each two-monomer propagator involves four endpoints, and thus, its diagram has four vertices which are connected to each other. The placement of the auxiliary variables in Eqs. (2.8) and (2.9) (one in the bra and one in the ket of each two-monomer factor, except for the end monomers which include one auxiliary and three endpoint variables) leads to the skewness observed in Fig. 1, i.e., the end units are represented by trapezoids, while all intermediate units appear as parallelograms. It is easy to check that the diagram of Eq. (2.8) is the mirror image of that corresponding to Eq. (2.9).

Next, consider the following Trotter factorization of the short-time evolution operator:

which is obtained from the left-to-right factorization given in Eq. (2.6) by swapping the monomers within each pair. This factorization leads to the following short-time propagator:

The diagram corresponding to this propagator is shown in Fig. 2. All units besides the first and last now have propagators connecting two endpoints in the bra (or ket) to two auxiliary variables in the ket (or bra). As a consequence, the skewness of the two-monomer propagators is mostly removed, i.e., all intermediate units are now represented by rectangles.

Each two-monomer unit is a system of M² states and thus is easily diagonalizable. To avoid introducing additional Trotter errors, each of the two-monomer propagators is to be obtained from the exact eigenstates of the corresponding Hamiltonian.

In this section, we present the structure of the MPI for the propagator at a time t = NΔt,

from the short-time propagators developed in Sec. II C. Inserting the resolution of identity N−1 times, the propagator becomes

Using the short time propagator according to the sequential factorization of Eq. (2.9), we obtain the following discretized path integral expression, which contains an auxiliary variable for each propagator of all monomers, except the first and last units:

All monomers in Eq. (3.3) (except for those at the two ends of the chain) involve 2N − 1 path integral variables. As a result, the MPI scheme corresponding to Eq. (3.3) requires the storage of M^2N−1-dimensional arrays, making the algorithm considerably more demanding compared to the MPI method with diagonal coupling operators, which requires the storage of M^N−1 paths for the same task. A diagrammatic representation of the path integral expression with this factorization is given in the left panel of Fig. 3. It is seen that each variable is connected to four other variables which link three adjacent monomers. We emphasize that different horizontal lengths of the parallelograms and trapezoids do not indicate different time lengths.

Next, consider the swapped pair factorization of Eq. (2.11), and its reverse order, which corresponds to the left/right reflection of the propagator diagram. Using an even value of N and alternating between these two propagators in adjacent path integral time steps, we obtain an MPI expression that corresponds to the diagram shown in the right panel of Fig. 3. Even though there are still 2N − 1 variables for each nonterminal monomer, one observes that the path integral variables corresponding to propagator endpoints (solid circles in Fig. 3) now have just three connections, which link each such variable to variables within the same two-monomer pair.

Furthermore, note that the first two scalene trapezoids (bottom right of Fig. 3), which correspond to the 0–1 and 1–2 propagators of the pair AB, form a larger isosceles trapezoid, and the vertical line connecting the adjacent vertices links the common endpoints of these propagators. Their contribution to the path integral expression has the form

Most importantly, the variables $n1A$ and $n1B$ do not appear in other factors in the path integral expression. Consequently, the sum over these variables in Eq. (3.4) is equal to the identity and may be removed, leading to the simpler form

which corresponds to the larger isosceles trapezoid in Fig. 3 that describes the double time step propagator for the AB pair. It is easy to see that all overlapping vertical edges of rectangles in the diagrammatic representation of Fig. 3 may also be removed this way, linking adjacent rectangles into squares that correspond to two-monomer propagators with time step 2Δt. This procedure allows elimination of all original path integral variables $nkα,k=1,…,N−1$⁠. The remaining sums in the path integral expression involve the N auxiliary variables of intermediate monomers and only ½N − 1 variables for the end units.

Renaming the auxiliary variables $ñk−1,k$ as $nkα,k=1,…,N$ for all nonterminal units gives rise to the new path integral expression,

The diagrammatic representation of Eq. (3.6), along with its symmetrization through reverse ordering, is shown in Fig. 4 for the case of five monomers with N = 4.

Equation (3.6) is in the desired form, where the number of path integral variables is almost the same as that in the case of coupling operators that are diagonal in a single-particle basis. We emphasize that the algorithm may be used in any single-monomer basis, e.g., either in the basis of single-particle eigenstates or in terms of localized states in which the single-monomer Hamiltonians are not diagonal.

In Sec. IV, we describe the implementation of the MPI algorithm to the present case of general, nondiagonal coupling operators, along with a factorization that dramatically reduces the computational effort in the sequential linking of monomer paths.

The path integral variables of each monomer in Eq. (3.6) are coupled only to those of adjacent monomers. Thus, this expression is suitable for the MPI decomposition, i.e., the sums over the variables for each monomer can be evaluated in a sequential manner.

The MPI calculation begins with the construction of an array (vector) $RA$ that contains the discrete paths of monomer A. Since the endpoints are fixed, there are $(MA)N$ such paths. Each element of this vector is to be linked to each path of monomer B via a matrix that contains the propagators for the AB pair. The process appears to scale as $(MB)N(MA)12N−1$⁠, but the path variables of monomer A can be summed one at a time. The scheme is an extension of a simpler procedure developed for MPI with diagonal couplings.⁴¹ 

The factorized MPI algorithm starts by summing over the first variable $n2A$ of monomer A. To eliminate this variable, it is necessary to construct an array of all variables to which $n2A$ couples, i.e., the first four variables of monomer B. Thus, the first step involves the multiplication,

Equation (4.1) involves M_A operations, and the resulting array has length $(MB)4MA$⁠. Next, the variable $n4A$ is eliminated after introducing the unit B variables to which it couples,

This procedure is repeated until the last variable of unit A is eliminated and leads to the array

Equations (4.1)–(4.3) are a total of ½N − 1 iterations, each of which involves a sum of M_A terms. The maximum storage is the array of all monomer B paths, which has $(MB)N$ elements. Thus, the total number of operations is $12N−1MA(MB)N$⁠.

The next step is to sum the paths of monomer B, after linking them to those of monomer C. We begin by eliminating the first variable by performing the multiplication,

We then eliminate variables $n2B$ and $n3B$ simultaneously,

and this process is repeated until all variables of monomer B have been dealt with. Each step in this process involves a sum of (M_B)² terms. Note that the required storage involves two arrays in each step. If M_B = M_C ≡ M, these arrays have length M^N.

The steps of Eqs. (4.4) and (4.5) are repeated until all the next-to-last monomer is reached, and that is linked to the terminal unit by analogous procedures. Summing all the elements of the resulting array produces the desired propagator amplitude.

We illustrate the methodology with application to two TLS chains (i.e., M_α = 2).

In the first model, the units are coupled spins. The single-particle Hamiltonians are given by

in the basis of the two eigenstates of σ_z for each unit, which are denoted as R and L (with the corresponding eigenvalues ±1), where $ℏΩ$ is half of the tunneling splitting. The coupling between adjacent TLSs has the form

with $cx=0.2ℏΩ$⁠, $cy=0.3ℏΩ$⁠, and $cz=0.4ℏΩ$⁠. Initially, all spins are in the R eigenstate of σ_z (with an eigenvalue +1). We present results with up to r = 31 spins.

Figure 5(a) shows the survival amplitude of the initial state,

for spin chains of length r = 3, 7, 25, and 31. Pronounced recurrences with interference patterns are seen in the shortest chain. The recurrences in the survival probability become progressively weaker as the chain length is increased.

Figure 5(b) shows the survival probability of the initial state, the probability of all spins tunneling to the L state,

and the probability of just the first spin tunneling, while the rest remain in their initial state,

The survival probability is seen to decay rapidly but shows recurrences of the decreasing magnitude. The probability of all spins flipping exhibits peaks which alternate with the peaks of the survival probability.

The second model focuses on the probability of energy transfer out of the excited state of the first unit in a chain of coupled chromophores. The monomer Hamiltonians are given by

where the energy gap $E1α−E0α≡2ℏΩ$⁠. The coupling has the Frenkel exciton form,

with $J=−0.47ℏΩ$⁠. These parameters are characteristic of conducting polymers.⁴² Initially, the first monomer is excited, while all other units are in the ground state.

Figure 6 shows monomer populations in two chains composed of 5 and 25 units. As time progresses, population transfer to successive units is observed. The time length over which the excitation resides on a particular unit exhibits a slow increase with the monomer distance from the initially populated edge unit. The shorter chain shows population accumulation on the last unit and a reflection of the excitation.

The MPI formulation allows fully quantum mechanical calculations in systems composed of segments of a one-dimensional topology and offers linear scaling with the system size. While this scaling is similar to that of DMRG methods, the unique advantage of the path integral formulation is its ability to yield zero- or finite-temperature properties within the same framework, and the possibility of including any number of molecular vibrations and the dissipative effects of phonons, which do not appear practical by other means. In this sense, the MPI decomposition promises to extend the capabilities of the fully quantum mechanical QuAPI methodology to systems composed of a large number of units, where the Hilbert space (which for a system of r monomers has dimension M^r) corresponding to the electronic states would be prohibitively large. The extension of the present scheme to the finite-temperature density matrices with explicit treatment of molecular vibrations will be described in future work.

Currently, the main limitation of the MPI methodology is the potentially large number of Feynman paths for each monomer, which grows exponentially with the number of path integral time steps. The calculations presented in Sec. V demonstrate that the converged results on state populations may be obtained over times that are sufficiently long to capture the important features of the process under investigation. Furthermore, a host of strategies that exploits the memory quenching role of dissipative interactions^{20,29–31,43} may be utilized in the MPI context to substantially decrease the number of paths, and other strategies may also be fruitful.

Another potentially very useful future development of the MPI methodology is to combine it with the QCPI formulation, which would allow the rigorous treatment of electronic-nuclear dynamics within each monomer, without restrictions to harmonic models of vibrational degrees of freedom. This work will be the subject of future work in our group.

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