Leonel
Varvelo
,
Jacob K.
Lynd
and
Doran I. G.
Bennett
*
Department of Chemistry, Southern Methodist University, PO Box 750314, Dallas, TX, USA. E-mail: doranb@smu.edu
First published on 31st May 2021
Excited state carriers, such as excitons, can diffuse on the 100 nm to micron length scale in molecular materials but only delocalize over short length scales due to coupling between electronic and vibrational degrees-of-freedom. Here, we leverage the locality of excitons to adaptively solve the hierarchy of pure states equations (HOPS). We demonstrate that our adaptive HOPS (adHOPS) methodology provides a formally exact and size-invariant (i.e., ) scaling algorithm for simulating mesoscale quantum dynamics. Finally, we provide proof-of-principle calculations for exciton diffusion on linear chains containing up to 1000 molecules.
Simulating exciton transport dynamics in heterogeneous materials on the 10 nm to 1 µm length scale (i.e., the mesoscale) remains an outstanding theoretical challenge. Organic semiconductors often combine close intermolecular packing with correspondingly large coupling between electronic states (V) on adjacent molecules and large intramolecular electron-vibrational coupling (λ).3 Perturbative equations-of-motion, such as Förster theory, can be convenient for simulating large aggregates, but are not applicable when V and λ are comparable in magnitude. Similarly, in the absence of a clear separation of timescales between vibrational and electronic degrees-of-freedom, Markovian equations-of-motion, such as Redfield theory,4 struggle to capture the rich dynamics of excitation transport. There are a variety of non-perturbative, non-Markovian equations-of-motion, such as multi-layer multi-configuration time-dependent Hartree (ML-MCTDH),5,6 time-evolving density operator with orthogonal polynomials (TEDOPA),7 hierarchically-coupled equations-of-motion (HEOM),8 and quasi-adiabatic path integrals (QUAPI).9 All of these techniques, however, share an exponential scaling of computational complexity with the number of molecules. While efficient and parallelized implementations of formally exact methods have been developed – for example, distributed memory HEOM10,11 – the exponential scaling severely limits even high-performance simulations of molecular aggregates.
Recently, there have been a few notable developments towards highly-scalable equations-of-motion for exciton dynamics. Modular path integrals12,13 provide a dramatic reduction in computational cost of QUAPI, but retain an overall linear scaling with the number of molecules and are most efficient when molecules exhibit only nearest-neighbor coupling. Dissipation-assisted matrix product factorization (DAMPF)14 extends TEDOPA to efficiently describe large numbers of vibrational degrees-of-freedom (>10) on each molecule, but it maintains between a quadratic and cubic scaling with the number of molecules. For both modular path integrals and DAMPF, the residual scaling makes it challenging to apply these methods to mesoscale calculations containing thousands to millions of molecules. Indeed, any density matrix approach will suffer from residual scaling with system size at long times due to the spread of ensemble population density across molecules.
Stochastic simulations, which decompose the ensemble into a collection of excited trajectories, can enable calculations on arbitrarily large molecular aggregates, even at long time. Delocalized kinetic Monte Carlo15 and the kinetic Monte Carlo version of generalized Förster theory16 are stochastic approaches that calculate the rate of transport between clusters of strongly interacting molecules and can be readily extended to mesoscale calculations. Both of these methods, however, use a perturbative approximation to partition state space and calculate rates between adjacent spatial regions. The development of a non-perturbative, non-Markovian approach for mesoscale simulations would provide an important benchmark for new equations-of-motion and could offer insight into processes with debated mechanisms, such as charge separation in organic photovoltaic materials.15,17–19
Here, we present a non-perturbative, non-Markovian, and arbitrarily scalable stochastic method for simulating exciton transport. First, we introduce the Hamiltonian considered and our base equation-of-motion, the hierarchy of pure states (HOPS).20 Next, we discuss locality in HOPS calculations and present an algorithm for constructing an adaptive basis. Finally, we present proof-of-concept calculations using the adaptive HOPS (adHOPS) equation-of-motion that demonstrate both its accuracy and size-invariant (i.e., ) scaling for large molecular aggregates.
![]() | (1) |
![]() | (2) |
![]() | (3) |
αn(t) = gne−γnt/ℏ | (4) |
![]() | (5) |
![]() | (6) |
The NMQSD equation is formally exact and is equivalent to solving Feynman path integrals with the Feynman–Vernon influence functional,22 but the functional derivative in the last term makes direct solution of the stochastic trajectories impractical except in special cases.
The hierarchy of pure states (HOPS) equations provide a numerically tractable version of NMQSD by rewriting the functional derivative as a set of coupled differential equations.20 Briefly, the sum of integrals over a functional derivative in the final term of the NMQSD equation is defined as a sum of first order auxiliary wave functions:
![]() | (7) |
![]() | (8) |
![]() | (9) |
We can improve convergence with ensemble size by using the non-linear HOPS equation, which describes the time-evolution of a normalizable stochastic wave function. We rewrite the reduced system density matrix (eqn (5)) in terms of a normalized wave function and the norm contribution
![]() | (10) |
![]() | (11) |
![]() | (12) |
![]() | (13) |
We note that the non-linear HOPS equation ensures that the contribution of each wave function is normalized in the reduced density matrix, but it time-evolves the non-normalized physical wave function. In the following, we will drop the explicit zn,t dependence from the wave function for simplicity .
The HOPS equations are a numerically convenient, formally exact expression for exciton dynamics in small molecular aggregates. Moreover, the calculations are ‘embarrassingly’ (also called ‘perfectly’) parallel24 due to the independence of individual trajectories, and, as a result, HOPS ensembles can be computed using thousands of CPUs simultaneously without loss of efficiency. The application of HOPS to large molecular aggregates, however, is limited by the scaling of the HOPS basis with the number of molecules. It is convenient to think of HOPS calculations as depending on two basis sets: the state basis and the auxiliary basis
. The complete state basis is a finite set of vectors that span the Hilbert space of the system, while the complete auxiliary basis is composed of an infinite set of auxiliary wave functions indexed by vectors
. To construct a finite auxiliary basis, the infinite hierarchy must be truncated. Here, we employ the common triangular truncation condition which limits the auxiliary basis to those wave functions with index vectors
that have a sum of elements less than a preselected bound kmax
. If we assume one independent environment per state, then the number of auxiliary wave functions included in the triangular truncation scales as
which gives an overall
scaling for large aggregates. While convergence as a function of kmax is guaranteed, the requisite number of auxiliary wave functions is often prohibitive.
![]() | (14) |
αdl,n(t) = (Re[gn] + isgn(t)Im[gn])e−γn|t|/ℏ | (15) |
The discontinuity in the correlation function introduces a numerically inconvenient infinitely high-frequency component to the stochastic noise trajectories zn,t.
We ameliorate this problem by redefining the positive-time correlation function in terms of two continuous exponential functions α(t) = α0,n(t) + αmark,n(t) where
α0,n(t) = gne−γn|t|/ℏ | (16) |
αmark,n(t) = −iIm[gn]e−γmark|t|/ℏ. | (17) |
The definition of αmark,n(t) ensures the imaginary component of the total correlation is 0 when t = 0, and it also provides for a smooth transition back to the naive correlation function α0,n(t) on a finite timescale given by ℏ/γmark. Except where otherwise noted, γmark = 500 cm−1 for all calculations presented in the main text, because this was sufficiently fast to ensure the Markovian timescale had no influence on the calculated dynamics.
Due to the extremely rapid timescale on which αmark,n(t) decays, this mode is Markovian, and high-lying contributions to the hierarchy can be neglected. In the following, we will only include the first order terms associated with these Markovian modes and neglect these terms in our discussion of the auxiliary wave functions forming the hierarchy. This can be viewed as a smoothing of the noise trajectories (zn,t) on timescales fast compared to all other dynamics.
This problem can be avoided entirely by using a different spectral density which more naturally accounts for the short-time imaginary component of the correlation function: for example, the recently reported alternative to the Drude–Lorentz oscillator with improved low-temperature behavior.26
Here, we develop an adaptive solution to the HOPS equation-of-motion that achieves size-invariant computational scaling (i.e., scaling) for calculations of large molecular aggregates. We first establish a normalized non-linear HOPS equation, which ensures that the magnitude of derivative terms does not diverge with increasing depth of the hierarchy. We then illustrate how locality appears within the hierarchy of auxiliary wave functions, with a particular emphasis on the connection between locality and the flux between neighboring auxiliary wave functions. Finally, we present an adaptive algorithm for the normalized non-linear HOPS equation that satisfies a user-selected bound on the absolute derivative error.
To enforce the normalization of the physical wave function, we rewrite the non-linear HOPS equation in terms of a normalized physical wave function. Starting with eqn (11), dividing all wave functions by the norm of the physical wave function, taking the derivative, and expanding terms gives
![]() | (18) |
![]() | (19) |
In the non-linear HOPS equation, the magnitude of the auxiliary wave functions grows with increasing auxiliary index. The basic HOPS terminator for a hierarchy with a single thermal environment
![]() | (20) |
![]() | (21) |
![]() | (22) |
Fig. 2 shows that in HOPS calculations localization in the physical wave function induces localization in the hierarchy. By ‘localization in the hierarchy,’ we specifically refer to clustering of amplitude in a small set of auxiliary wave functions in a way that depends on the position of the excitation in the physical wave function. In Fig. 2a, an excitation begins on the middle site of a five pigment chain and then jumps between site 3 and site 2. Fig. 2b plots the norm of auxiliary wave functions associated with site 2 and site 3; each plot is labeled by an auxiliary vector index . For example, [0,1,0,0,0] (Fig. 2b, first column) is the first order auxiliary wave function associated with the thermal environment on the second site. The occupation of the auxiliary wave functions (black lines, Fig. 2b) track with the population of the physical wave function – i.e., the auxiliary wave functions associated with the thermal environment of site 2 (first column Fig. 2b) are only occupied when site 2 is occupied in the physical wave function (shaded region).
The locality in the HOPS hierarchy can be understood in terms of the balance of flux terms in the normalized non-linear HOPS equation. First, every auxiliary wave function is damped (first line of eqn (21), ) and, therefore, has zero amplitude without a continuous source term. The fundamental source term is the physical wave function which is the lowest order of the hierarchy. The flux of amplitude towards higher-lying auxiliary wave functions arises from the second line of eqn (21)
. For system-bath coupling operators that are site projection operators
, the flux towards higher-lying auxiliaries only arises when there is amplitude on the associated site of the lower auxiliary wave function. Moreover, the auxiliary wave functions are localized by the same dynamics that localize the physical wave function. As a result, the localized auxiliary wave functions only contribute amplitude to higher-lying auxiliary wave functions with an index that differs by
in a site (n) with non-zero amplitude. Thus, the locality of the physical wave function results in preferential population of specific auxiliary wave functions.
Our adHOPS algorithm neither assumes nor imposes locality. Rather, the adaptive basis takes advantage of whatever locality arises during the simulation. If the full hierarchy is required to satisfy the derivative error bound, adHOPS smoothly reverts to a HOPS calculation. As a result, adHOPS remains formally exact – the adaptive basis represents a time-dependent truncation of hierarchy elements, and δ, like kmax, is a convergence parameter.
We note that the current adHOPS algorithm makes use of two approximations: first, the spectral density is assumed to be over-damped (e.g., Drude–Lorentz), which allows for a consistent normalization of the hierarchy elements. Second, the system-bath coupling operator is assumed to be a site-projection operator (n = |n〉〈n|); in other words, we assume that each molecule has an independent vibrational environment.
One persistent challenge for numerical implementations of formally exact methods is demonstrating the calculations are converged to the exact answer. In hierarchical methods, calculations must be converged with respect to the auxiliary basis which is defined in the triangular truncation condition by the maximum hierarchy level considered (kmax). In HOPS, the criterion for convergence is that kmaxγ/ℏ ≫ ωs, where ωs is the characteristic frequency of the system.20 Because the full auxiliary basis scales as , it is often impractical to systematically check convergence for sufficiently large values of kmax. Though our adHOPS method was inspired by localization, we find that it naturally incorporates a dynamic filtering scheme that dramatically improves the scaling of the auxiliary basis with kmax even when the exciton is fully delocalized. Fig. 4a compares the full (black) and adaptive (green) HOPS dynamics with increasing coupling (V). By V = 5λ the oscillations in the site 3 population report a wave function that is coherently oscillating across 5 sites. Fig. 4b shows the corresponding size of the auxiliary basis as a function of kmax. In all cases the adaptive auxiliary basis (green line) increases much more slowly than the full auxiliary basis (black line).
Another perpetual challenge for formally exact methods is their intractable computational scaling with the number of molecules. In HOPS calculations this arises from the scaling of the auxiliary basis. Fig. 5a compares the full (black line) size of the state (top) and auxiliary (bottom) basis to the average size of the adaptive basis (colored lines) as a function of the number of molecules in a linear chain. Compared to the full auxiliary basis, the size of the adaptive auxiliary basis scales favorably with respect to the number of molecules in the aggregate. Moreover, both the auxiliary and state bases in adHOPS calculations show a plateau beyond a threshold size of the linear chain (Npig > N*), indicating the onset of size invariant scaling. In the ESI,† we compare the CPU time required for full and adaptive HOPS calculations (V = 50 cm−1). We find that adaptive calculations are faster than full calculations starting around Npig = 10, and we also demonstrate the onset of size invariance (i.e., ) scaling of CPU time for large aggregates. In other words, increasing the number of pigments beyond a threshold size does not increase the computational expense of an adHOPS calculation. Thus, for localized excitons, the size invariance of adHOPS allows for calculations on scales that were previously unachievable for formally exact methods.
Our adaptive HOPS algorithm offers a computationally tractable approach for formally exact calculations of mesoscale quantum dynamics. As a proof-of-concept, we demonstrate the ability to simulate exciton diffusion on a linear chain of 103 molecules within the formally exact framework of adHOPS (Fig. 5b). Exciton diffusion is a common experimental observable extracted from non-linear microscopies2 but is challenging to simulate on long length scales.35–37 Using adHOPS, simulating exciton diffusion in a linear chain of 103 pigments is computationally tractable because for V = 100 cm−1 it requires, on average, less than 2 × 103 auxiliary wave functions and 20 pigment states. The corresponding HOPS simulation would require an auxiliary basis containing more than 1023 auxiliary wave functions.
1. Is a formally exact solution to the time-evolution of a quantum state coupled to a non-Markovian thermal reservoir,
2. Is embarrassingly (or ‘perfectly’) parallel,24 and
3. Achieves size-invariant (i.e., ) scaling for molecular aggregates that are substantially larger than the exciton delocalization extent in the material.
This combination of properties allows us to perform non-perturbative, non-Markovian simulations involving an arbitrary number of pigments in physically relevant parameter regimes, thus laying the foundation for mesoscale quantum dynamics simulations of excited-state carriers in molecular materials. Currently, our adaptive algorithm assumes that each pigment has an independent thermal environment composed of overdamped vibrations, but future developments will allow for a broader class of mechanisms involving high-frequency intra-molecular vibrations38–40 and Peierls-type electron-vibration coupling.3 Looking forward, we believe that adHOPS provides a promising new direction for simulations of a broad range of organic semiconductors including photosynthetic membranes,16,41 molecular thin films,42,43 and organic photovoltaic heterojunctions.17–19,44
Footnote |
† Electronic supplementary information (ESI) available: PDF contains the explanation of the statistical, hierarchy, and time-step errors; the expressions for evaluating derivative error terms; description of the adaptive error distribution in the ensemble; CPU timing data for adHOPS calculations. See DOI: 10.1039/d1sc01448j |
This journal is © The Royal Society of Chemistry 2021 |